home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdsample.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-28  |  7.8 KB  |  248 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdsample.c    2.20  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.  * grdsample reads a grdfile and evaluates the grid at new grid positions
  9.  * specified by new dx/dy values using a 2-D Taylor expansion of order 3.
  10.  * In order to evaluate derivatives along the edges of the surface, I assume 
  11.  * natural bi-cubic spline conditions, i.e. both the second and third normal 
  12.  * derivatives are zero, and that the dxdy derivative in the corners are zero, too.
  13.  *
  14.  * Author:    Paul Wessel
  15.  * Date:    19-JUL-1989
  16.  * Revised:    6-JAN-1990    PW: Updated to v.2.0
  17.  */
  18.  
  19. #include "gmt.h"
  20. #include "gmt_bcr.h"
  21.  
  22. float *a, *b;
  23.  
  24. main (argc, argv)
  25. int argc;
  26. char **argv; {
  27.     int i, j, ij, one, pad[4];
  28.     
  29.     BOOLEAN error = FALSE, greenwich = FALSE, offset = FALSE, bilinear = FALSE;
  30.     BOOLEAN area_set = FALSE, n_set = FALSE, inc_set = FALSE, geo = FALSE, toggle = FALSE, global_grid;
  31.     
  32.     double *lon, lat, dx2, dy2;
  33.     
  34.     char *infile, *outfile;
  35.     
  36.     struct GRD_HEADER grd_a, grd_b;
  37.  
  38.     argc = gmt_begin (argc, argv);
  39.     
  40.     infile = outfile = CNULL;
  41.     
  42.     grd_init (&grd_b, argc, argv, FALSE);
  43.  
  44.     for (i = 1; i < argc; i++) {
  45.         if (argv[i][0] == '-') {
  46.             switch (argv[i][1]) {
  47.         
  48.                 /* Common parameters */
  49.             
  50.                 case 'R':
  51.                 case 'V':
  52.                 case '\0':
  53.                     error += get_common_args (argv[i], &grd_b.x_min, &grd_b.x_max, &grd_b.y_min, &grd_b.y_max);
  54.                     break;
  55.                 
  56.                 /* Supplemental parameters */
  57.             
  58.                 case 'F':
  59.                     offset = TRUE;
  60.                     break;
  61.                 case 'G':
  62.                     outfile = &argv[i][2];
  63.                     break;
  64.                 case 'N':
  65.                     sscanf (&argv[i][2], "%d/%d", &grd_b.nx, &grd_b.ny);
  66.                     if (grd_b.ny == 0) grd_b.ny = grd_b.nx;
  67.                     n_set = TRUE;
  68.                     break;
  69.                 case 'I':
  70.                     gmt_getinc (&argv[i][2], &grd_b.x_inc, &grd_b.y_inc);
  71.                     inc_set = TRUE;
  72.                     break;
  73.                 case 'L':
  74.                     geo = TRUE;
  75.                     break;
  76.                 case 'Q':
  77.                     bilinear = TRUE;
  78.                     break;
  79.                 case 'T':    /* Convert from pixel file <-> gridfile */
  80.                     toggle = TRUE;
  81.                     break;
  82.                 default:
  83.                     error = TRUE;
  84.                     gmt_default_error (argv[i][1]);
  85.                     break;
  86.             }
  87.         }
  88.         else 
  89.             infile = argv[i];
  90.     }
  91.     
  92.     if (argc == 1 || gmt_quick) {
  93.         fprintf(stderr, "grdsample %s - Resample a gridded file onto a new grid\n\n", GMT_VERSION);
  94.         fprintf(stderr, "usage: grdsample <old_grdfile> -G<new_grdfile> [-F] [-I<dx>[m|c][/<dy>[m|c]]]\n");
  95.         fprintf(stderr, "    [-L] [-N<nx/ny>] [-Q] [-R<west/east/south/north>] [-T] [-V]\n");
  96.         
  97.         if (gmt_quick) exit (-1);
  98.         
  99.         fprintf(stderr, "    <old_grdfile> is data set to be resampled\n");
  100.         fprintf(stderr, "    -G sets the name of the interpolated output grdfile\n");
  101.         fprintf(stderr, "\n\tOPTIONS:\n");
  102.         fprintf(stderr, "    -F force pixel node registration  [Default is centered]\n");
  103.         fprintf(stderr, "    -I sets the grid spacing (dx, dy) for the new grid\n");
  104.         fprintf(stderr, "    -L means these are geographical grids\n");
  105.         fprintf(stderr, "    -N specifies number of columns (nx) and rows (ny) of new grid\n");
  106.         fprintf(stderr, "    -Q do quick, bilinear interpolation [Default is bicubic]\n");
  107.         fprintf(stderr, "    -R specifies a subregion [Default is old region]\n");
  108.         fprintf(stderr, "    -T Toggles between grid registration and pixel registration\n");
  109.         explain_option ('V');
  110.         fprintf(stderr, "       One only of -N -I must be specified\n");
  111.         
  112.         exit (-1);
  113.     }
  114.     
  115.     if (!infile) {
  116.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  117.         error++;
  118.     }
  119.     if (!outfile) {
  120.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", gmt_program);
  121.         error++;
  122.     }
  123.     if (!toggle) {
  124.         if (grd_b.x_min != grd_b.x_max && grd_b.y_min != grd_b.y_max) area_set = TRUE;
  125.         if (inc_set && n_set) {
  126.             fprintf (stderr, "%s: GMT SYNTAX ERROR:  Only one of -I, -N may be specified\n", gmt_program);
  127.             error++;
  128.         }
  129.         if (n_set && (grd_b.nx <= 0 || grd_b.ny <= 0)) {
  130.             fprintf (stderr, "%s: GMT SYNTAX ERROR -N:  Must specify positive integers\n", gmt_program);
  131.             error++;
  132.         }
  133.         if (inc_set && (grd_b.x_inc <= 0.0 || grd_b.y_inc <= 0.0)) {
  134.             fprintf (stderr, "%s: GMT SYNTAX ERROR -I:  Must specify positive increments\n", gmt_program);
  135.             error++;
  136.         }
  137.         if (!(inc_set || n_set)) {
  138.             fprintf (stderr, "%s: GMT SYNTAX ERROR:  One of -I, -N must be specified\n", gmt_program);
  139.             error++;
  140.         }
  141.     }
  142.     
  143.     if (error) exit (-1);
  144.     
  145.     if (read_grd_info (infile, &grd_a)) {
  146.         fprintf (stderr, "grdsample: Error opening file %s\n", infile);
  147.         exit (-1);
  148.     }
  149.     
  150.     if (toggle) {
  151.         offset = !grd_a.node_offset;    /* Change to the opposite of what it is */
  152.         grd_b.nx = (offset) ? grd_a.nx - 1 : grd_a.nx + 1;
  153.         grd_b.ny = (offset) ? grd_a.ny - 1 : grd_a.ny + 1;
  154.         area_set = inc_set = FALSE;
  155.     }
  156.     
  157.     a = (float *) memory (CNULL, (int)((grd_a.nx + 4) * (grd_a.ny + 4)), sizeof(float), "grdsample");
  158.     
  159.     global_grid = (fabs (grd_a.x_max - grd_a.x_min) == 360.0);
  160.  
  161.     if (area_set) {
  162.         if (grd_b.y_min < grd_a.y_min || grd_b.y_max > grd_a.y_max) {
  163.             fprintf (stderr, "grdsample:  Selected region exceeds the boundaries of the grdfile!\n");
  164.             exit (-1);
  165.         }
  166.         else if (!global_grid && (grd_b.x_min < grd_a.x_min || grd_b.x_max > grd_a.x_max)) {
  167.             fprintf (stderr, "grdsample:  Selected region exceeds the boundaries of the grdfile!\n");
  168.             exit (-1);
  169.         }
  170.     }
  171.     
  172.     if (!offset && !toggle) offset = grd_a.node_offset;
  173.     one = !offset;
  174.     grd_b.node_offset = offset;
  175.     
  176.     if (!area_set) {
  177.         grd_b.x_min = grd_a.x_min;
  178.         grd_b.x_max = grd_a.x_max;
  179.         grd_b.y_min = grd_a.y_min;
  180.         grd_b.y_max = grd_a.y_max;
  181.     }
  182.     
  183.     if (geo && grd_b.x_min < 0.0 && grd_b.x_max > 0.0)
  184.         greenwich = TRUE;
  185.     else if (geo && grd_b.x_max > 360.0) {
  186.         greenwich = TRUE;
  187.         grd_b.x_min -= 360.0;
  188.         grd_b.x_max -= 360.0;
  189.     }
  190.     if (inc_set) {
  191.         grd_b.nx = rint ((grd_b.x_max - grd_b.x_min) / grd_b.x_inc) + one;
  192.         grd_b.ny = rint ((grd_b.y_max - grd_b.y_min) / grd_b.y_inc) + one;
  193.         grd_b.x_inc = (grd_b.x_max - grd_b.x_min) / (grd_b.nx - one);
  194.         grd_b.y_inc = (grd_b.y_max - grd_b.y_min) / (grd_b.ny - one);
  195.     }
  196.     else {
  197.         grd_b.x_inc = (grd_b.x_max - grd_b.x_min) / (grd_b.nx - one);
  198.         grd_b.y_inc = (grd_b.y_max - grd_b.y_min) / (grd_b.ny - one);
  199.     }
  200.     
  201.     b = (float *) memory (CNULL, (int)(grd_b.nx * grd_b.ny), sizeof(float), "grdsample");
  202.     
  203.     if (gmtdefs.verbose) fprintf (stderr, "grdsample: New grid (%lg/%lg/%lg/%lg) nx = %d ny = %d dx = %lg dy = %lg\n",
  204.         grd_b.x_min, grd_b.x_max, grd_b.y_min, grd_b.y_max, grd_b.nx, grd_b.ny, grd_b.x_inc, grd_b.y_inc);
  205.     
  206.     pad[0] = pad[1] = pad[2] = pad[3] = 2;    /* Leave room for 2 empty boundary rows/cols */
  207.     
  208.     if (read_grd (infile, &grd_a, a, grd_a.x_min, grd_a.x_max, grd_a.y_min, grd_a.y_max, pad, FALSE)) {
  209.         fprintf (stderr, "grdsample: Error reading file %s\n", infile);
  210.         exit (-1);
  211.     }
  212.  
  213.     /* Initialize bcr structure:  */
  214.  
  215.     bcr_init (&grd_a, pad, bilinear);
  216.  
  217.     /* Set boundary conditions  */
  218.     
  219.     set_bcr_boundaries (&grd_a, a);
  220.     
  221.     /* Precalculate longitudes */
  222.     
  223.     dx2 = 0.5 * grd_a.x_inc;
  224.     dy2 = 0.5 * grd_a.y_inc;
  225.     lon = (double *) memory (CNULL, grd_b.nx, sizeof (double), "grdsample");
  226.     for (i = 0; i < grd_b.nx; i++) {
  227.         lon[i] = grd_b.x_min + (i * grd_b.x_inc) + ((offset) ? dx2 : 0.0);
  228.         if (geo && greenwich && lon[i] > 180.0) lon[i] -= 360.0;
  229.     }
  230.  
  231.     for (j = ij = 0; j < grd_b.ny; j++) {
  232.         lat = grd_b.y_max - (j * grd_b.y_inc);
  233.         if (offset) lat -= dy2;
  234.         for (i = 0; i < grd_b.nx; i++, ij++) b[ij] = get_bcr_z (&grd_a, lon[i], lat, a);
  235.     }
  236.     
  237.     if (write_grd (outfile, &grd_b, b, 0.0, 0.0, 0.0, 0.0, pad, FALSE)) {
  238.         fprintf (stderr, "grdsample: Error writing file %s\n", outfile);
  239.         exit (-1);
  240.     }
  241.     
  242.     free ((char *)a);
  243.     free ((char *)b);
  244.     free ((char *)lon);
  245.     
  246.     gmt_end (argc, argv);
  247. }
  248.