home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdinfo.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-25  |  5.6 KB  |  188 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdinfo.c    2.22  5/25/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.  * grdinfo reads one or more grd file and [optionally] prints out various
  9.  * statistics like mean/standard deviation and median/scale
  10.  *
  11.  * Author:    Paul Wessel
  12.  * Date:    4-JAN-1991-1995
  13.  * Version:    3.0    Based on 2.x
  14.  */
  15.  
  16. #include "gmt.h"
  17.  
  18. float *a;
  19. char *type[2] = { "Normal", "Pixel"};
  20.  
  21. main (argc, argv)
  22. int argc;
  23. char **argv; {
  24.     int nfiles = 0, k, i, j, i_min, i_max, nm, n_nan = 0, n, dummy[4];
  25.     
  26.     BOOLEAN error = FALSE, l1 = FALSE, l2 = FALSE, quick = TRUE, find_max = FALSE;
  27.     
  28.     double x_min, y_min, z_min, x_max, y_max, z_max;
  29.     double mean, median, sum2, stdev, scale, rms, x;
  30.     
  31.     char file[BUFSIZ];
  32.     
  33.     struct GRD_HEADER grd;
  34.     
  35.     argc = gmt_begin (argc, argv);
  36.  
  37.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  38.  
  39.     for (i = 1; i < argc; i++) {
  40.         if (argv[i][0] == '-') {
  41.             switch (argv[i][1]) {
  42.                 /* Common parameters */
  43.             
  44.                 case '\0':
  45.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  46.                     break;
  47.                 
  48.                 /* Supplemental parameters */
  49.             
  50.                 case 'M':
  51.                     quick = FALSE;
  52.                     find_max = TRUE;
  53.                     break;
  54.                 case 'L':
  55.                     quick = FALSE;
  56.                     if (argv[i][2] == 0 || argv[i][2] == '2')
  57.                         l2 = TRUE;
  58.                     else if (argv[i][2] == '1')
  59.                         l1 = TRUE;
  60.                     else {
  61.                         error = TRUE;
  62.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -L option:  Choose between -L1 or -L2\n", gmt_program);
  63.                     }
  64.                     break;
  65.                 default:
  66.                     error = TRUE;
  67.                     gmt_default_error (argv[i][1]);
  68.                     break;
  69.             }
  70.         }
  71.         else
  72.             nfiles ++;
  73.     }
  74.     
  75.     if (argc == 1 || gmt_quick) {
  76.         fprintf (stderr, "grdinfo %s - Extract information from netCDF grdfiles\n\n", GMT_VERSION);
  77.         fprintf (stderr, "usage: grdinfo <grdfiles> [-L1] [-L[2]] [-M]\n");
  78.         
  79.         if (gmt_quick) exit (-1);
  80.         
  81.         fprintf (stderr, "    <grdfiles> may be one or more netCDF grdfiles\n");
  82.         fprintf (stderr, "\n\tOPTIONS:\n");
  83.         fprintf (stderr, "    -L1 will report median/L1-scale of data set\n");
  84.         fprintf (stderr, "    -L[2] will report mean/standard deviation/rms of data set\n");
  85.         fprintf (stderr, "    -M will search for location of min/max z-values\n");
  86.         exit (-1);
  87.     }
  88.     
  89.     if (nfiles == 0) {
  90.         fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one or more input files\n", gmt_program);
  91.         error++;
  92.     }
  93.     if (error) exit (-1);
  94.  
  95.     a = (float *) memory (CNULL, 1, sizeof (float), "grdinfo");
  96.  
  97.     for (k = 1; k < argc; k++) {    /* Loop over arguments, skip options */
  98.     
  99.         if (argv[k][0] == '-') continue;
  100.         
  101.         strcpy (file, argv[k]);
  102.         for (j = 0; file[j]; j++) if (file[j] == '=') file[j] = 0;
  103.         if (strcmp (file, "=") && access (file, R_OK)) {
  104.             fprintf (stderr, "grdinfo:  File %s not found\n", file);
  105.             continue;
  106.         }
  107.  
  108.         grd_init (&grd, argc, argv, FALSE);
  109.  
  110.         if (read_grd_info (argv[k], &grd)) {
  111.             fprintf (stderr, "grdinfo: Error opening file %s\n", file);
  112.             continue;
  113.         }
  114.         if (grd.z_min == grd.z_max && grd_i_format) quick = FALSE, find_max = TRUE;
  115.         
  116.         printf ("%s: Title: %s\n", file, grd.title);
  117.         printf ("%s: Command: %s\n", file, grd.command);
  118.         printf ("%s: Remark: %s\n", file, grd.remark);
  119.         printf ("%s: %s node registation used\n", file, type[grd.node_offset]);
  120.         printf ("%s: grdfile format # %d\n", file, grd_i_format);
  121.         printf ("%s: x_min: %lg x_max: %lg x_inc: %lg units: %s nx: %d\n", file, grd.x_min, grd.x_max, grd.x_inc, grd.x_units, grd.nx);
  122.         printf ("%s: y_min: %lg y_max: %lg y_inc: %lg units: %s ny: %d\n", file, grd.y_min, grd.y_max, grd.y_inc, grd.y_units, grd.ny);
  123.         if (!quick) {
  124.             nm = grd.nx * grd.ny;
  125.             a = (float *) memory ((char *)a, nm, sizeof (float), "grdinfo");
  126.             if (read_grd (argv[k], &grd, a, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) continue;
  127.         
  128.             z_min = MAXDOUBLE;    z_max = -MAXDOUBLE;
  129.             mean = median = sum2 = 0.0;
  130.             i_min = i_max = 0;
  131.             n_nan = 0;
  132.             for (i = 0; i < nm; i++) {
  133.                 if (bad_float ((double)a[i])) {
  134.                     n_nan++;
  135.                     continue;
  136.                 }
  137.                 if (find_max) {
  138.                     if (a[i] < z_min) {
  139.                         z_min = a[i];
  140.                         i_min = i;
  141.                     }
  142.                     if (a[i] > z_max) {
  143.                         z_max = a[i];
  144.                         i_max = i;
  145.                     }
  146.                 }
  147.                 if (l2) {
  148.                     mean += a[i];
  149.                     sum2 += a[i]*a[i];
  150.                 }
  151.             }
  152.         
  153.             x_min = grd.x_min + (i_min % grd.nx) * grd.x_inc;
  154.             y_min = grd.y_max - (i_min / grd.nx) * grd.y_inc;
  155.             x_max = grd.x_min + (i_max % grd.nx) * grd.x_inc;
  156.             y_max = grd.y_max - (i_max / grd.nx) * grd.y_inc;
  157.         }
  158.         
  159.         if (find_max) 
  160.             printf ("%s: zmin: %lg (at %lg/%lg) zmax: %lg (at %lg/%lg) units: %s\n", file, z_min, x_min, y_min, z_max, x_max, y_max, grd.z_units);
  161.         else
  162.             printf ("%s: zmin: %lg zmax: %lg units: %s\n", file, grd.z_min, grd.z_max, grd.z_units);
  163.  
  164.         printf ("%s: scale_factor: %lg add_offset: %lg\n", file, grd.z_scale_factor, grd.z_add_offset);
  165.         if (n_nan) printf ("%s: %d nodes set to NaN\n", file, n_nan);
  166.         if (l1) {
  167.             qsort ((char *)a, nm, sizeof (float), comp_float_asc);
  168.             n = nm - n_nan;
  169.             median = (n%2) ? a[n/2] : 0.5*(a[n/2-1] + a[n/2]);
  170.             for (i = 0; i < n; i++) a[i] = fabs (a[i] - median);
  171.             qsort ((char *)a, n, sizeof (float), comp_float_asc);
  172.             scale = (n%2) ? 1.4826 * a[n/2] : 0.7413 * (a[n/2-1] + a[n/2]);
  173.             printf ("%s: median: %lg scale: %lg\n", file, median, scale);
  174.         }
  175.         if (l2) {
  176.             x = nm - n_nan;
  177.             stdev = sqrt((x*sum2 - mean*mean)/(x*(x-1)));
  178.             rms = sqrt (sum2 / x);
  179.             mean /= x;
  180.             printf ("%s: mean: %lg stdev: %lg rms: %lg\n", file, mean, stdev, rms);
  181.         }
  182.     }
  183.     
  184.     free ((char *)a);
  185.     
  186.     gmt_end (argc, argv);
  187. }
  188.