home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdlandmask.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-07  |  12.1 KB  |  336 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdlandmask.c    3.14  07 Aug 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.  * grdlandmask defines a grid based on region and xinc/yinc values,
  9.  * reads a shoreline data base, 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:    P. Wessel
  14.  * Date:    23-Sep-1994
  15.  * Version:    3.0
  16.  */
  17.  
  18. #include "gmt.h"
  19. #include "gmt_shore.h"        /* Defines shore structures */
  20.  
  21. char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};
  22.  
  23. struct SHORE c;
  24.  
  25. float *data;
  26.  
  27. main (argc, argv)
  28. int argc;
  29. char **argv;
  30. {
  31.  
  32.     int    i, j, k, ij, bin, ind, nm, np, side, i_min, i_max, j_min, j_max, nx1, ny1;
  33.     int    one_or_zero, base = 3, direction, inside = 1, max_level = MAX_LEVEL, dummy[4];
  34.     
  35.     BOOLEAN    error = FALSE, pixel = FALSE, dry_wet_only = FALSE, greenwich = FALSE, got_fh[2];
  36.     BOOLEAN temp_shift = FALSE;
  37.     
  38.     double    *x, *y, xmin, xmax, ymin, ymax, out_edge_in[5];
  39.     double xinc2, yinc2, min_area = 0.0, i_dx, i_dy, edge = 0.0, del_off;
  40.     
  41.     char *maskfile = CNULL, res = 'l', line[80], *ptr;
  42.     
  43.     struct GRD_HEADER header;
  44.     struct POL *p;
  45.     
  46.     argc = gmt_begin (argc, argv);
  47.         
  48.     grd_init (&header, argc, argv, FALSE);
  49.     
  50.     memset ((char *)out_edge_in, 0, 5 * sizeof (double));    /* Default "wet" value = 0 */
  51.     out_edge_in[1] = out_edge_in[3] = 1.0;            /* Default "dry" value = 1 */
  52.  
  53.     /* Check command line arguments */
  54.     
  55.     for (i = 1; i < argc; i++) {
  56.         if (argv[i][0] == '-') {
  57.             switch (argv[i][1]) {
  58.             
  59.                 /* Common parameters */
  60.             
  61.                 case 'R':
  62.                 case 'V':
  63.                 case ':':
  64.                 case '\0':
  65.                     error += get_common_args (argv[i], &header.x_min, &header.x_max, &header.y_min, &header.y_max);
  66.                     break;
  67.                 
  68.                 /* Supplemental parameters */
  69.                 
  70.                 case 'A':
  71.                     j = sscanf (&argv[i][2], "%lf/%d", &min_area, &max_level);
  72.                     if (j == 1) max_level = MAX_LEVEL;
  73.                     break;
  74.                 case 'D':
  75.                     res = argv[i][2];
  76.                     switch (res) {
  77.                         case 'f':
  78.                             base = 0;
  79.                             break;
  80.                         case 'h':
  81.                             base = 1;
  82.                             break;
  83.                         case 'i':
  84.                             base = 2;
  85.                             break;
  86.                         case 'l':
  87.                             base = 3;
  88.                             break;
  89.                         case 'c':
  90.                             base = 4;
  91.                             break;
  92.                         default:
  93.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Unknown modifier %c\n", gmt_program, argv[i][2]);
  94.                             break;
  95.                     }
  96.                     break;
  97.                 case 'N':
  98.                     strcpy (line, &argv[i][2]);
  99.                     if (line[strlen(line)-1] == 'o') { /* Edge is considered outside */
  100.                         inside = 2;
  101.                         line[strlen(line)-1] = 0;
  102.                     }
  103.                     ptr = strtok (line, "/");
  104.                     j = 0;
  105.                     while (j < 5 && ptr) {
  106.                         out_edge_in[j] = (ptr[0] == 'N' || ptr[0] == 'n') ? gmt_NaN : atof (ptr);
  107.                         ptr = strtok (CNULL, "/");
  108.                         j++;
  109.                     }
  110.                     if (!(j == 2 || j == 5)) {
  111.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Specify 2 or 5 arguments\n", gmt_program);
  112.                         exit (-1);
  113.                     }
  114.                     dry_wet_only = (j == 2);
  115.                     break;
  116.                 case 'F':
  117.                     pixel = TRUE;
  118.                     break;
  119.                 case 'G':
  120.                     maskfile = &argv[i][2];
  121.                     break;
  122.                 case 'I':
  123.                     gmt_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
  124.                     break;
  125.                 default:
  126.                     error = TRUE;
  127.                     gmt_default_error (argv[i][1]);
  128.                     break;
  129.             }
  130.         }
  131.     }
  132.  
  133.         find_resolutions (got_fh);      /* Check to see if full &| high resolution coastlines are installed */
  134.  
  135.     if (argc == 1 || gmt_quick) {
  136.         fprintf (stderr, "grdlandmask %s - Create \"wet-dry\" mask grdfile from shoreline data base\n\n", GMT_VERSION);
  137.         fprintf (stderr, "usage: grdlandmask -G<mask_grd_file> -I<xinc[m|c]>[/<yinc>[m|c]] -R<west/east/south/north>\n");
  138.         fprintf (stderr, "\t[-A<min_area>[/<max_level>]] [-D<resolution>] [-F] [-N<maskvalues>[o]] [-V]\n\n");
  139.         
  140.         if (gmt_quick) exit (-1);
  141.         
  142.         fprintf (stderr, "\t-G Specify file name for output mask grd file.\n");
  143.         fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
  144.         fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to indicate minutes [or seconds],\n");
  145.         explain_option ('R');
  146.         fprintf (stderr, "\n\tOPTIONS:\n");
  147.         fprintf (stderr, "\t-A features smaller than <min_area> (in km^2) or of higher level (0-4) than <max_level>\n");
  148.         fprintf (stderr, "\t   will be skipped [0/4] (see pscoast for details)]\n");
  149.         fprintf (stderr, "    -D Choose one of the following resolutions:\n");
  150.         if (got_fh[0]) fprintf (stderr, "       f - full resolution (may be very slow for large regions)\n");
  151.         if (got_fh[1]) fprintf (stderr, "       h - high resolution (may be slow for large regions)\n");
  152.         fprintf (stderr, "       i - intermediate resolution\n");
  153.         fprintf (stderr, "       l - low resolution [Default]\n");
  154.         fprintf (stderr, "       c - crude resolution, for tasks that need crude continent outlines only\n");
  155.         fprintf (stderr, "\t-F Force pixel registration for output grid [Default is gridline orientation]\n");
  156.         fprintf (stderr, "\t-N gives values to use if a node is outside or inside a feature.\n");
  157.         fprintf (stderr, "\t   Append o to let feature boundary be considered outside [Default is inside].\n");
  158.         fprintf (stderr, "\t   Specify this information using 1 of 2 formats:\n");
  159.         fprintf (stderr, "\t   -N<wet>/<dry>.\n");
  160.         fprintf (stderr, "\t   -N<ocean>/<land>/<lake>/<island>/<pond>.\n");
  161.         fprintf (stderr, "\t   NaN is a valid entry.  Default values are 0/1/0/1/0 (i.e., 0/1)\n");
  162.         explain_option ('V');
  163.         exit(-1);
  164.     }
  165.  
  166.     if (!project_info.region_supplied) {
  167.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  168.         error++;
  169.     }
  170.     if (header.x_inc <= 0.0 || header.y_inc <= 0.0) {
  171.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", gmt_program);
  172.         error = TRUE;
  173.     }
  174.     if (!maskfile) {
  175.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", gmt_program);
  176.         error = TRUE;
  177.     }
  178.     
  179.     if (error) exit (-1);
  180.  
  181.     if (header.x_min < 0.0 && header.x_max < 0.0) {    /* Shift longitudes */
  182.         temp_shift = TRUE;
  183.         header.x_min += 360.0;
  184.         header.x_max += 360.0;
  185.     }
  186.     
  187.     if (dry_wet_only) {
  188.         out_edge_in[3] = out_edge_in[5] = out_edge_in[7] = out_edge_in[1];
  189.         out_edge_in[4] = out_edge_in[8] = out_edge_in[0];
  190.         out_edge_in[6] = out_edge_in[2];
  191.     }
  192.  
  193.     if (gmt_init_shore (res, &c, header.x_min, header.x_max, header.y_min, header.y_max))  {
  194.         fprintf (stderr, "grdlandmask: %s resolution shoreline data base not installed\n", shore_resolution[base]);
  195.         exit (-1);
  196.     }
  197.     
  198.     if (gmtdefs.verbose && dry_wet_only) {
  199.         fprintf (stderr, "grdlandmask: Nodes in water will be set to ");
  200.         (bad_float (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[0]);
  201.         fprintf (stderr, "grdlandmask: Nodes on land will be set to ");
  202.         (bad_float (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[1]);
  203.     }
  204.     else if (gmtdefs.verbose && !dry_wet_only) {
  205.         fprintf (stderr, "grdlandmask: Nodes in the oceans will be set to ");
  206.         (bad_float (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[0]);
  207.         fprintf (stderr, "grdlandmask: Nodes on land will be set to ");
  208.         (bad_float (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[1]);
  209.         fprintf (stderr, "grdlandmask: Nodes in lakes will be set to ");
  210.         (bad_float (out_edge_in[2])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[2]);
  211.         fprintf (stderr, "grdlandmask: Nodes in islands will be set to ");
  212.         (bad_float (out_edge_in[3])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[3]);
  213.         fprintf (stderr, "grdlandmask: Nodes in ponds will be set to ");
  214.         (bad_float (out_edge_in[4])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", out_edge_in[4]);
  215.     }
  216.     
  217.     i_dx = 1.0 / header.x_inc;
  218.     i_dy = 1.0 / header.y_inc;
  219.     one_or_zero = (pixel) ? 0 : 1;
  220.     del_off = (pixel) ? 0.0 : 0.5;
  221.     header.nx = rint((header.x_max - header.x_min) * i_dx) + one_or_zero;
  222.     header.ny = rint((header.y_max - header.y_min) * i_dy) + one_or_zero;
  223.     header.node_offset = pixel;
  224.     xinc2 = (header.node_offset) ? 0.5 * header.x_inc : 0.0;
  225.     yinc2 = (header.node_offset) ? 0.5 * header.y_inc : 0.0;
  226.     nm = header.nx * header.ny;
  227.     data = (float *) memory (CNULL, nm, sizeof(float), "grdlandmask");
  228.     x = (double *) memory (CNULL, header.nx, sizeof(double), "grdlandmask");
  229.     y = (double *) memory (CNULL, header.ny, sizeof(double), "grdlandmask");
  230.     for (i = 0; i < header.nx; i++) x[i] = header.x_min + i * header.x_inc + xinc2;
  231.     for (j = 0; j < header.ny; j++) y[j] = header.y_max - j * header.y_inc - yinc2;
  232.     nx1 = header.nx - 1;    ny1 = header.ny - 1;
  233.  
  234.     if (header.x_min < 0.0 && header.x_max > 0.0) {    /* Must shift longitudes */
  235.         greenwich = TRUE;
  236.         edge = header.x_min;
  237.     }
  238.     
  239.     map_getproject ("x1d");    /* Fake linear projection */
  240.     map_setup (header.x_min, header.x_max, header.y_min, header.y_max);
  241.     
  242.     for (ind = 0; ind < c.nb; ind++) {    /* Loop over necessary bins only */
  243.         
  244.         bin = c.bins[ind];
  245.         if (gmtdefs.verbose) fprintf (stderr, "grdlandmask: Working on block # %5d\r", bin);
  246.  
  247.         gmt_get_shore_bin (ind, &c, min_area, max_level);
  248.         
  249.         if (c.ns == 0) { /* Initialize all nodes to lowest node level */
  250.         
  251.             k = MIN (MIN (c.node_level[0], c.node_level[1]) , MIN (c.node_level[2], c.node_level[3]));
  252.             i_min = MAX (0, floor (fmod (c.lon_sw - header.x_min + 360.0, 360.0) * i_dx) + del_off);
  253.             if (i_min > nx1) i_min = 0;
  254.             i_max = MIN (nx1, ceil (fmod (c.lon_sw + c.bsize - header.x_min + 360.0, 360.0) * i_dx) + del_off);
  255.             if (i_max <= 0) i_max = nx1;
  256.             j_min = MAX (0, floor ((header.y_max - c.lat_sw - c.bsize) * i_dy) - del_off);
  257.             j_max = MIN (ny1, ceil ((header.y_max - c.lat_sw) * i_dy) - del_off);
  258.             for (j = j_min; j <= j_max; j++) for (i = i_min, ij = j * header.nx + i_min; i <= i_max; i++, ij++)
  259.                 data[ij] = k;
  260.         
  261.             continue;    /* No polygon, done with bin */
  262.         }
  263.  
  264.         /* Must use polygons.  Go in both directions to cover both land and sea */
  265.         
  266.         for (direction = -1; direction < 2; direction += 2) {
  267.         
  268.             np = gmt_assemble_shore (&c, direction, TRUE, greenwich, edge, &p);
  269.             
  270.             for (k = 0; k < np; k++) {
  271.         
  272.                 /* Find min/max of polygon */
  273.             
  274.                 if (greenwich && p[k].lon[0] < edge) p[k].lon[0] += 360.0;
  275.                 xmin = xmax = p[k].lon[0];
  276.                 ymin = ymax = p[k].lat[0];
  277.                 for (i = 1; i < p[k].n; i++) {
  278.                     if (greenwich && p[k].lon[i] < edge) p[k].lon[i] += 360.0;
  279.                     if (p[k].lon[i] < xmin) xmin = p[k].lon[i];
  280.                     if (p[k].lon[i] > xmax) xmax = p[k].lon[i];
  281.                     if (p[k].lat[i] < ymin) ymin = p[k].lat[i];
  282.                     if (p[k].lat[i] > ymax) ymax = p[k].lat[i];
  283.                 }
  284.                 i_min = MAX (0, floor (fmod (xmin - header.x_min + 360.0, 360.0) * i_dx) + del_off);
  285.                 if (i_min > nx1) i_min = 0;
  286.                 i_max = MIN (nx1, ceil (fmod (xmax - header.x_min + 360.0, 360.0) * i_dx) + del_off);
  287.                 if (i_max <= 0) i_max = nx1;
  288.                 j_min = MAX (0, floor ((header.y_max - ymax) * i_dy) - del_off);
  289.                 j_max = MIN (ny1, ceil ((header.y_max - ymin) * i_dy) - del_off);
  290.             
  291.                 for (j = j_min; j <= j_max; j++) {
  292.                     for (i = i_min; i <= i_max; i++) {
  293.                     
  294.                         if ((side = non_zero_winding (x[i], y[j], p[k].lon, p[k].lat, p[k].n)) < inside) continue;    /* Outside */
  295.                 
  296.                         /* Here, point is inside, we must assign value */
  297.                 
  298.                         ij = j * header.nx + i;
  299.                         if (p[k].level > data[ij]) data[ij] = p[k].level;
  300.                     }
  301.                 }
  302.             }
  303.         
  304.             free_polygons (p, np);
  305.             free ((char *)p);
  306.         }
  307.         
  308.         free_shore (&c);
  309.     }
  310.     
  311.     shore_cleanup (&c);
  312.     
  313.     for (ij = 0; ij < nm; ij++) {
  314.         k = rint (data[ij]);
  315.         data[ij] = out_edge_in[k];
  316.     }
  317.     
  318.     if (temp_shift) {
  319.         header.x_min -= 360.0;
  320.         header.x_max -= 360.0;
  321.     }
  322.  
  323.     if (write_grd (maskfile, &header, data, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  324.         fprintf (stderr, "grdlandmask: Error writing file %s\n", maskfile);
  325.         exit (-1);
  326.     }
  327.  
  328.     if (gmtdefs.verbose) fprintf (stderr, "\ngrdlandmask: Done!\n");
  329.     
  330.     free ((char *)data);
  331.     free ((char *)x);
  332.     free ((char *)y);
  333.     
  334.     gmt_end (argc, argv);
  335. }
  336.