home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdproject.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-30  |  8.2 KB  |  250 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdproject.c    2.19  30 Jun 1995
  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.  * grdproject reads a geographical grdfile and evaluates the grid at new grid positions
  9.  * specified by the map projection and new dx/dy values using a weighted average of all
  10.  * points within the search radius. Optionally, grdproject may perform the inverse
  11.  * transformation, going from rectangular coordinates to geographical.
  12.  *
  13.  * Author:    Paul Wessel
  14.  * Date:    6-JAN-1991-1995
  15.  * Ver:        2.0
  16.  *
  17.  */
  18.  
  19. #include "gmt.h"
  20.  
  21. float *geo, *rect, *weight;
  22.  
  23. main (argc, argv)
  24. int argc;
  25. char **argv; {
  26.     int i, dpi, nx, ny, dummy[4];
  27.     
  28.     BOOLEAN error = FALSE, inverse = FALSE, n_set = FALSE, set_n = FALSE;
  29.     BOOLEAN d_set = FALSE, e_set = FALSE, map_center = FALSE, offset = FALSE;
  30.     
  31.     double w, e, s, n, x_inc = 0.0, y_inc = 0.0, max_radius = 0.0;
  32.     double xmin, xmax, ymin, ymax;
  33.     
  34.     char *infile, *outfile;
  35.     
  36.     struct GRD_HEADER g_head, r_head;
  37.     
  38.     argc = gmt_begin (argc, argv);
  39.     
  40.     infile = outfile = CNULL;
  41.     w = e = s = n = 0.0;
  42.     dummy[3] = dummy[2] = dummy[1] = dummy[0] = 0;
  43.     
  44.     for (i = 1; i < argc; i++) {
  45.         if (argv[i][0] == '-') {
  46.             switch (argv[i][1]) {
  47.                 /* Common parameters */
  48.             
  49.                 case 'J':
  50.                 case 'R':
  51.                 case 'V':
  52.                 case '\0':
  53.                     error += get_common_args (argv[i], &w, &e, &s, &n);
  54.                     break;
  55.                 
  56.                 /* Supplemental parameters */
  57.                 
  58.                 case 'D':
  59.                     gmt_getinc (&argv[i][2], &x_inc, &y_inc);
  60.                     d_set = TRUE;
  61.                     break;
  62.                 case 'E':
  63.                     dpi = atoi (&argv[i][2]);
  64.                     e_set = TRUE;
  65.                     break;
  66.                 case 'F':
  67.                     offset = TRUE;
  68.                     break;
  69.                 case 'G':
  70.                     outfile = &argv[i][2];
  71.                     break;
  72.                 case 'I':
  73.                     inverse = TRUE;
  74.                     break;
  75.                 case 'M':
  76.                     map_center = TRUE;
  77.                     if (argv[i][2] == 'm') gmtdefs.measure_unit = 2;    /* Select meters */
  78.                     break;
  79.                 case 'N':
  80.                     sscanf (&argv[i][2], "%d/%d", &nx, &ny);
  81.                     if (ny == 0) ny = nx;
  82.                     n_set = TRUE;
  83.                     break;
  84.                 case 'S':
  85.                     max_radius = atof (&argv[i][2]);
  86.                     break;
  87.                 default:
  88.                     error = TRUE;
  89.                     gmt_default_error (argv[i][1]);
  90.                     break;
  91.             }
  92.         }
  93.         else 
  94.             infile = argv[i];
  95.     }
  96.     
  97.     if ((d_set + e_set + n_set) == 0) n_set = set_n = TRUE;
  98.  
  99.     if (argc == 1 || gmt_quick) {
  100.         fprintf (stderr, "grdproject %s - Project geographical grid to/from rectangular grid\n\n", GMT_VERSION);
  101.         fprintf(stderr, "usage: grdproject <in_grdfile> -J<parameters> -R<west/east/south/north>\n");
  102.         fprintf(stderr, "    [-D<dx[m|c]>[/dy[m|c]]] [-E<dpi>] [-F] [-G<out_grdfile>] [-I] [-M]\n");
  103.         fprintf(stderr, "    [-N<nx/ny>] [-S<radius>] [-V]\n\n");
  104.         
  105.         if (gmt_quick) exit (-1);
  106.         
  107.         fprintf(stderr, "    <in_grdfile> is data set to be transformed\n");
  108.         explain_option ('J');
  109.         explain_option ('R');
  110.         fprintf(stderr, "\n\tOPTIONS:\n");
  111.         fprintf(stderr, "    -D sets the grid spacing for the new grid\n");
  112.         fprintf(stderr, "    -E sets dpi for output grid\n");
  113.         fprintf(stderr, "    -F force pixel node registration  [Default is centered]\n");
  114.         fprintf(stderr, "    -G name of output grid\n");
  115.         fprintf(stderr, "    -I Inverse transformation from rectangular to geographical\n");
  116.         fprintf (stderr, "    -M coordinates relative to projection center [Default is relative to lower left corner]\n");
  117.         fprintf(stderr, "    -N sets the number of nodes for the new grid\n");
  118.         fprintf(stderr, "       Only one of -D, -E, and -N can be specified!\n");
  119.         fprintf(stderr, "       If none are specified, nx,ny of the input file is used\n");
  120.         fprintf(stderr, "    -S sets the search radius in projected units [Default avoids aliasing]\n");
  121.         explain_option ('V');
  122.         
  123.         exit (-1);
  124.     }
  125.     
  126.     if (!infile) {
  127.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  128.         error++;
  129.     }
  130.     if (!outfile) {
  131.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", gmt_program);
  132.         error++;
  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 (max_radius < 0.0) {
  139.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S:  Max search radius, if specified, must be positive\n", gmt_program);
  140.         error++;
  141.     }
  142.     if ((d_set + e_set + n_set) != 1) {
  143.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify only one of -D, -E, or -N\n", gmt_program);
  144.         error++;
  145.     }
  146.     if (d_set && (x_inc <= 0.0 || y_inc <= 0.0)) {
  147.         fprintf (stderr, "%s: GMT SYNTAX ERROR -D option.  Must specify positive increment(s)\n", gmt_program);
  148.         error++;
  149.     }
  150.     if (n_set && !set_n && (nx <= 0 || ny <= 0)) {
  151.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  Must specify positive integers\n", gmt_program);
  152.         error++;
  153.     }
  154.     if (e_set && dpi <= 0) {
  155.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option.  Must specify positive dpi\n", gmt_program);
  156.         error++;
  157.     }
  158.     
  159.     if (error) exit (-1);
  160.  
  161.     map_setup (w, e, s, n);
  162.     
  163.     xmin = (map_center) ? project_info.xmin - project_info.x0 : project_info.xmin;
  164.     xmax = (map_center) ? project_info.xmax - project_info.x0 : project_info.xmax;
  165.     ymin = (map_center) ? project_info.ymin - project_info.y0 : project_info.ymin;
  166.     ymax = (map_center) ? project_info.ymax - project_info.y0 : project_info.ymax;
  167.     
  168.     grd_init (&r_head, argc, argv, FALSE);
  169.     grd_init (&g_head, argc, argv, FALSE);
  170.  
  171.     if (inverse) {    /* Transforming from rectangular projection to geographical */
  172.     
  173.         if (!project_info.region) d_swap (s, e);  /* Got w/s/e/n, make into w/e/s/n */
  174.  
  175.         g_head.x_min = w;    g_head.x_max = e;    g_head.y_min = s;    g_head.y_max = n;
  176.     
  177.         if (read_grd_info (infile, &r_head)) {
  178.             fprintf (stderr, "grdproject: Error opening file %s\n", infile);
  179.             exit (-1);
  180.         }
  181.         rect = (float *) memory (CNULL, (int)(r_head.nx * r_head.ny), sizeof (float), "grdproject");
  182.         if (read_grd (infile, &r_head, rect, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  183.             fprintf (stderr, "grdproject: Error reading file %s\n", infile);
  184.             exit (-1);
  185.         }
  186.     
  187.         if (set_n) {
  188.             nx = r_head.nx;
  189.             ny = r_head.ny;
  190.         }
  191.         grdproject_init (&g_head, x_inc, y_inc, nx, ny, dpi, offset);
  192.         geo = (float *) memory (CNULL, (int)(g_head.nx * g_head.ny), sizeof (float), "grdproject");
  193.         if (gmtdefs.verbose)
  194.             fprintf (stderr, "grdproject:  Transform (%lg/%lg/%lg/%lg) <-- (%lg/%lg/%lg/%lg)\n",
  195.                 g_head.x_min, g_head.x_max, g_head.y_min, g_head.y_max, xmin, xmax, ymin, ymax);
  196.  
  197.         grd_inverse (geo, &g_head, rect, &r_head, max_radius, map_center);
  198.             
  199.         if (write_grd (outfile, &g_head, geo, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  200.             fprintf (stderr, "grdproject: Error writing file %s\n", outfile);
  201.             exit (-1);
  202.         }
  203.         
  204.     }
  205.     else {    /* Forward projection from geographical to rectangular grid */
  206.     
  207.         if (read_grd_info (infile, &g_head)) {
  208.             fprintf (stderr, "grdproject: Error opening file %s\n", infile);
  209.             exit (-1);
  210.         }
  211.  
  212.         if (project_info.projection == MOLLWEIDE && (fabs (g_head.x_min - g_head.x_max) < 360.0 || (g_head.y_max - g_head.y_min) < 180.0)) {
  213.             fprintf (stderr, "grdproject: Selected projection works on global datasets only\n");
  214.             exit (-1);
  215.         }
  216.     
  217.         geo = (float *) memory (CNULL, (int)(g_head.nx * g_head.ny), sizeof (float), "grdproject");
  218.         if (read_grd (infile, &g_head, geo, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  219.             fprintf (stderr, "grdproject: Error reading file %s\n", infile);
  220.             exit (-1);
  221.         }
  222.     
  223.         r_head.x_min = xmin;    r_head.x_max = xmax;
  224.         r_head.y_min = ymin;    r_head.y_max = ymax;
  225.         if (set_n) {
  226.             nx = g_head.nx;
  227.             ny = g_head.ny;
  228.         }
  229.         
  230.         if (gmtdefs.verbose)
  231.             fprintf (stderr, "grdproject:  Transform (%lg/%lg/%lg/%lg) --> (%lg/%lg/%lg/%lg)\n",
  232.                 g_head.x_min, g_head.x_max, g_head.y_min, g_head.y_max, xmin, xmax, ymin, ymax);
  233.  
  234.         grdproject_init (&r_head, x_inc, y_inc, nx, ny, dpi, offset);
  235.         rect = (float *) memory (CNULL, (int)(r_head.nx * r_head.ny), sizeof (float), "grdproject");
  236.         grd_forward (geo, &g_head, rect, &r_head, max_radius, map_center);
  237.         
  238.         if (write_grd (outfile, &r_head, rect, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  239.             fprintf (stderr, "grdproject: Error writing file %s\n", outfile);
  240.             exit (-1);
  241.         }
  242.         
  243.     }
  244.     
  245.     free ((char *)geo);
  246.     free ((char *)rect);
  247.     
  248.     gmt_end (argc, argv);
  249. }
  250.