home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdmask.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-28  |  8.1 KB  |  254 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdmask.c    2.21  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.  * grdmask defines a grid based on region and xinc/yinc values,
  9.  * reads xy polygon files, and sets the grid nodes inside, on the
  10.  * boundary, and outside of the polygons to the user-defined values
  11.  * <in>, <on>, and <out>.  These may be any number, including NaN.
  12.  *
  13.  * Author:    Walter. H. F. Smith
  14.  * Date:    23-May-1991-1995
  15.  * Version:    2.0
  16.  */
  17.  
  18. #include "gmt.h"
  19.  
  20. float *data;
  21.  
  22. main (argc, argv)
  23. int argc;
  24. char **argv;
  25. {
  26.  
  27.     int    i, j, ij, nm, n_path, ix, iy, side, fno, n_files = 0, n_alloc = GMT_CHUNK, n_args;
  28.     int    one_or_zero, dummy[4];
  29.     
  30.     BOOLEAN    error = FALSE, done, nofile = TRUE, multi_segments = FALSE, pixel = FALSE;
  31.     
  32.     double    *x, *y, xx, yy, xmin, xmax, ymin, ymax, xy[2], out_edge_in[3];
  33.     double xinc2, yinc2;
  34.     
  35.     char *maskfile = CNULL, line[512], *ptr, *more, EOL_flag = '>';
  36.     
  37.     FILE *fp;
  38.     
  39.     struct GRD_HEADER header;
  40.     
  41.     argc = gmt_begin (argc, argv);
  42.         
  43.     grd_init (&header, argc, argv, FALSE);
  44.     
  45.     out_edge_in[2] = 1.0;            /* Default inside value */
  46.     out_edge_in[0] = out_edge_in[1] = 0.0;    /* Default outside and edge value */
  47.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  48.  
  49.     /* Check command line arguments */
  50.     
  51.     for (i = 1; i < argc; i++) {
  52.         if (argv[i][0] == '-') {
  53.             switch (argv[i][1]) {
  54.             
  55.                 /* Common parameters */
  56.             
  57.                 case 'R':
  58.                 case 'V':
  59.                 case ':':
  60.                 case '\0':
  61.                     error += get_common_args (argv[i], &header.x_min, &header.x_max, &header.y_min, &header.y_max);
  62.                     break;
  63.                 
  64.                 /* Supplemental parameters */
  65.                 
  66.                 case 'N':
  67.                     strcpy (line, &argv[i][2]);
  68.                     ptr = strtok (line, "/");
  69.                     for (j = 0; j < 3; j++) {
  70.                         out_edge_in[j] = (ptr[0] == 'N' || ptr[0] == 'n') ? gmt_NaN : atof (ptr);
  71.                         ptr = strtok (CNULL, "/");
  72.                     }
  73.                     break;
  74.                 case 'F':
  75.                     pixel = TRUE;
  76.                     break;
  77.                 case 'G':
  78.                     maskfile = &argv[i][2];
  79.                     break;
  80.                 case 'I':
  81.                     gmt_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
  82.                     break;
  83.                 case 'M':
  84.                     multi_segments = TRUE;
  85.                     if (argv[i][2]) EOL_flag = argv[i][2];
  86.                     break;
  87.                 default:
  88.                     error = TRUE;
  89.                     gmt_default_error (argv[i][1]);
  90.                     break;
  91.             }
  92.         }
  93.         else
  94.             n_files++;
  95.     }
  96.  
  97.     if (argc == 1 || gmt_quick) {
  98.         fprintf (stderr, "grdmask %s - Create mask grdfile from polygons\n\n", GMT_VERSION);
  99.         fprintf (stderr, "usage: grdmask [<pathfiles>] -G<mask_grd_file> -I<xinc[m|c]>[/<yinc>[m|c]]\n");
  100.         fprintf (stderr, "\t-R<west/east/south/north> [-F] [-N<out>/<edge>/<in>] [-M[<flag>]] [-V] [-:]\n\n");
  101.         
  102.         if (gmt_quick) exit (-1);
  103.         
  104.         fprintf (stderr, "\tpathfiles is one or more polygon files\n");
  105.         fprintf (stderr, "\t-G Specify file name for output mask grd file.\n");
  106.         fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
  107.         fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to indicate minutes [or seconds].\n");
  108.         explain_option ('R');
  109.         fprintf (stderr, "\n\tOPTIONS:\n");
  110.         fprintf (stderr, "\t-F Force pixel registration for output grid [Default is gridline orientation]\n");
  111.         fprintf (stderr, "\t-M Input files each consist of multiple segments separated by one record\n");
  112.         fprintf (stderr, "\t    whose first character is <flag> [>].\n");
  113.         fprintf (stderr, "\t-N sets values to use if point is outside, on the path, or inside.\n");
  114.         fprintf (stderr, "\t   NaN is a valid entry.  Default values are 0/0/1.\n");
  115.         explain_option ('V');
  116.         explain_option (':');
  117.         exit(-1);
  118.     }
  119.  
  120.     if (!project_info.region_supplied) {
  121.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  122.         error++;
  123.     }
  124.     if (header.x_inc <= 0.0 || header.y_inc <= 0.0) {
  125.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", gmt_program);
  126.         error = TRUE;
  127.     }
  128.     if (!maskfile) {
  129.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", gmt_program);
  130.         error = TRUE;
  131.     }
  132.     
  133.     if (error) exit (-1);
  134.  
  135.     one_or_zero = (pixel) ? 0 : 1;
  136.     header.nx = rint((header.x_max - header.x_min) / header.x_inc) + one_or_zero;
  137.     header.ny = rint((header.y_max - header.y_min) / header.y_inc) + one_or_zero;
  138.     header.node_offset = pixel;
  139.     xinc2 = 0.5 * header.x_inc;
  140.     yinc2 = 0.5 * header.y_inc;
  141.     nm = header.nx * header.ny;
  142.     data = (float *) memory (CNULL, nm, sizeof(float), "grdmask");
  143.     x = (double *) memory (CNULL, n_alloc, sizeof(double), "grdmask");
  144.     y = (double *) memory (CNULL, n_alloc, sizeof(double), "grdmask");
  145.  
  146.     if (gmtdefs.verbose) {
  147.         fprintf (stderr, "grdmask: Nodes completely outside the polygons will be set to ");
  148.         (bad_float (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[0]);
  149.         fprintf (stderr, "grdmask: Nodes completely inside the polygons will be set to ");
  150.         (bad_float (out_edge_in[2])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[2]);
  151.         fprintf (stderr, "grdmask: Nodes on the polygons boundary will be set to ");
  152.         (bad_float (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[1]);
  153.     }
  154.     
  155.     /* Initialize all nodes to the 'outside' value */
  156.     
  157.     for (ij = 0; ij < nm; ij++) data[ij] = out_edge_in[0];
  158.     
  159.     ix = (gmtdefs.xy_toggle) ? 1 : 0;    iy = 1 - ix;        /* Set up which columns have x and y */
  160.  
  161.     if (n_files > 0)
  162.         nofile = FALSE;
  163.     else
  164.         n_files = 1;
  165.     n_args = (argc > 1) ? argc : 2;
  166.     
  167.     done = FALSE;
  168.     
  169.     for (fno = 1; !done && fno < n_args; fno++) {    /* Loop over all input files */
  170.         if (!nofile && argv[fno][0] == '-') continue;
  171.         if (nofile) {
  172.             fp = stdin;
  173.             done = TRUE;
  174.         }
  175.         else if ((fp = fopen (argv[fno], "r")) == NULL) {
  176.             fprintf (stderr, "grdmask: Cannot open file %s\n", argv[fno]);
  177.             continue;
  178.         }
  179.  
  180.         if (!nofile && gmtdefs.verbose) fprintf (stderr, "grdmask: Working on file %s\n", argv[fno]);
  181.         
  182.         if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fp);
  183.         
  184.         more = fgets (line, 512, fp);
  185.         while (more) {
  186.             n_path = 0;
  187.             xmin = ymin = 1.0e+100;
  188.             xmax = ymax = -1.0e+100;
  189.     
  190.             while ((more && !multi_segments) || (more && multi_segments && line[0] != EOL_flag)) {
  191.         
  192.                 sscanf(line, "%lf %lf", &xy[ix], &xy[iy]);
  193.                 x[n_path] = xy[0];
  194.                 y[n_path] = xy[1];
  195.     
  196.                 if (x[n_path] > xmax) xmax = x[n_path];
  197.                 if (x[n_path] < xmin) xmin = x[n_path];
  198.                 if (y[n_path] > ymax) ymax = y[n_path];
  199.                 if (y[n_path] < ymin) ymin = y[n_path];
  200.                 n_path++;
  201.         
  202.                 if (n_path == (n_alloc-1)) {    /* n_alloc-1 since we may need 1 more to close polygon */
  203.                     n_alloc += GMT_CHUNK;
  204.                     x = (double *) memory ((char *)x, n_alloc, sizeof(double), "grdmask");
  205.                     y = (double *) memory ((char *)y, n_alloc, sizeof(double), "grdmask");
  206.                 }
  207.                 more = fgets (line, 512, fp);
  208.             }
  209.     
  210.             if (x[n_path-1] != x[0] || y[n_path-1] != y[0]) {
  211.                 x[n_path] = x[0];
  212.                 y[n_path] = y[0];
  213.                 n_path++;
  214.             }
  215.  
  216.             for (j = 0; j < header.ny; j++) {
  217.         
  218.                 yy = header.y_max - j * header.y_inc;
  219.                 if (pixel) yy -= yinc2;
  220.                 if (yy < ymin || yy > ymax) continue;    /* Point outside, no need to assign value */
  221.             
  222.                 for (i = 0; i < header.nx; i++) {
  223.                     xx = header.x_min + i * header.x_inc;
  224.                     if (pixel) xx += xinc2;
  225.                     if (xx < xmin || xx > xmax) continue;    /* Point outside, no need to assign value */
  226.                 
  227.                     if ((side = non_zero_winding(xx, yy, x, y, n_path)) == 0) continue;    /* Outside */
  228.                 
  229.                     /* Here, point is inside or on edge, we must assign value */
  230.                 
  231.                     data[j*header.nx+i] = out_edge_in[side];
  232.                 }
  233. #ifdef DEBUG
  234.                 fprintf (stderr, "j = %d\r", j);
  235. #endif
  236.             }
  237.             if (more) more = fgets (line, 512, fp);
  238.         }
  239.         
  240.         if (fp != stdin) fclose (fp);
  241.     }
  242.     
  243.     if (write_grd (maskfile, &header, data, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  244.         fprintf (stderr, "grdmask: Error writing file %s\n", maskfile);
  245.         exit (-1);
  246.     }
  247.  
  248.     free ((char *)data);
  249.     free ((char *)x);
  250.     free ((char *)y);
  251.     
  252.     gmt_end (argc, argv);
  253. }
  254.