home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdcut.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  4.2 KB  |  141 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdcut.c    2.14  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.  * grdcut.c reads a grdfile and writes a portion within it
  9.  * to a new file.
  10.  *
  11.  * Author:    Walter Smith
  12.  * Date:    5 august, 1988
  13.  *
  14.  */
  15.  
  16. #include "gmt.h"
  17.  
  18. float    *grd;
  19.  
  20. main (argc, argv)
  21. int argc;
  22. char **argv;
  23. {
  24.     int    error = FALSE, global = FALSE;
  25.     int    i, nx_old, ny_old, nx_new, ny_new, one_or_zero, dummy[4];
  26.     double    w_new = 0.0, e_new = 0.0, s_new = 0.0, n_new = 0.0;
  27.     double    w_old, e_old, s_old, n_old;
  28.     char *grd_in, *grd_out;
  29.     struct GRD_HEADER header;
  30.  
  31.     argc = gmt_begin (argc, argv);
  32.     
  33.     grd_in = grd_out = CNULL;
  34.     
  35.     /* Check and interpret the command line arguments */
  36.     
  37.     for (i =1; i < argc; i++) {
  38.         if (argv[i][0] == '-') {
  39.             switch(argv[i][1]) {
  40.                 /* Common parameters */
  41.             
  42.                 case 'R':
  43.                 case 'V':
  44.                 case '\0':
  45.                     error += get_common_args (argv[i], &w_new, &e_new, &s_new, &n_new);
  46.                     break;
  47.                     
  48.                  case 'G':
  49.                      grd_out = &argv[i][2];
  50.                     break;
  51.                 default:        /* Options not recognized */
  52.                     error = TRUE;
  53.                     gmt_default_error (argv[i][1]);
  54.                     break;
  55.             }
  56.         }
  57.         else
  58.              grd_in = argv[i];
  59.     }
  60.     
  61.     if (argc == 1 || gmt_quick) {
  62.         fprintf (stderr,"grdcut %s - Extract subsets from grdfiles\n\n", GMT_VERSION);
  63.         fprintf (stderr, "usage: grdcut <input_grd> -G<output_grd> -R<west/east/south/north> [-V]\n");
  64.         
  65.         if (gmt_quick) exit (-1);
  66.         
  67.         fprintf (stderr, "    <input_grd> is file to extract a subset from.\n");
  68.         fprintf (stderr, "    -G specifies output grdfile\n");
  69.         explain_option ('R');
  70.         fprintf (stderr, "       Obviously, the WESN you specify must be within the WESN of the input file.\n");
  71.         fprintf (stderr, "       If in doubt, run grdinfo first and check range of old file.\n");
  72.         fprintf (stderr, "\n\tOPTIONS:\n");
  73.         explain_option ('V');
  74.         exit (-1);
  75.     }
  76.  
  77.     /* Check that the options selected make sense */
  78.     
  79.     if (!project_info.region_supplied) {
  80.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  81.         error++;
  82.     }
  83.     if (!grd_out) {
  84.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", gmt_program);
  85.         error++;
  86.     }
  87.     if (!grd_in) {
  88.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  89.         error++;
  90.     }
  91.     if (error) exit (-1);
  92.  
  93.     if (read_grd_info (grd_in, &header)) {
  94.         fprintf (stderr, "grdcut: Error opening file %s\n", grd_in);
  95.         exit (-1);
  96.     }
  97.         
  98.     if (s_new < header.y_min || s_new > header.y_max) error = TRUE;
  99.     if (n_new < header.y_min || n_new > header.y_max) error = TRUE;
  100.     global = (fabs (header.x_max - header.x_min) == 360.0);
  101.         
  102.     if ( !global && ((w_new < header.x_min) || (e_new > header.x_max)) ) error = TRUE;
  103.     if (error) {
  104.         fprintf (stderr, "grdcut: Subset exceeds data domain!\n");
  105.         exit (-1);
  106.     }
  107.  
  108.     grd_init (&header, argc, argv, TRUE);
  109.  
  110.     w_old = header.x_min;    e_old = header.x_max;
  111.     s_old = header.y_min;    n_old = header.y_max;
  112.     nx_old = header.nx;    ny_old = header.ny;
  113.     one_or_zero = (header.node_offset) ? 0 : 1;
  114.     nx_new = rint ((e_new - w_new) / header.x_inc) + one_or_zero;
  115.     ny_new = rint ((n_new - s_new) / header.y_inc) + one_or_zero;
  116.             
  117.     grd = (float *) memory (CNULL, nx_new * ny_new, sizeof (float), "grdcut");
  118.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  119.     if (read_grd (grd_in, &header, grd, w_new, e_new, s_new, n_new, dummy, FALSE)) {
  120.         fprintf (stderr, "grdcut: Error reading file %s\n", grd_in);
  121.         exit (-1);
  122.     }
  123.  
  124.     if (gmtdefs.verbose) {
  125.         fprintf (stderr, "File spec:\tW E S N dx dy nx ny:\n");
  126.         fprintf (stderr, "Old:\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%d\t%d\n",
  127.             w_old, e_old, s_old, n_old, header.x_inc, header.y_inc, nx_old, ny_old);
  128.         fprintf (stderr, "New:\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%.8lg\t%d\t%d\n",
  129.             w_new, e_new, s_new, n_new, header.x_inc, header.y_inc, nx_new, ny_new);
  130.     }
  131.  
  132.     if (write_grd (grd_out, &header, grd, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  133.         fprintf (stderr, "grdcut: Error writing file %s\n", grd_out);
  134.         exit (-1);
  135.     }
  136.     
  137.     free ( (char *) grd);
  138.     
  139.     gmt_end (argc, argv);
  140. }
  141.