home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / gmt_shore.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-11-20  |  27.5 KB  |  792 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)gmt_shore.c    3.9  26 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. #include "gmt.h"
  9. #include "gmt_shore.h"
  10.  
  11. int find_resolutions (fh)
  12. BOOLEAN fh[2];
  13. {
  14.     /* Checks if full and high resolution coastline files are available */
  15.  
  16.     char file[80];
  17. /*----------------------------------------------------------------------
  18.  *      Modifications made to permit the use of an environmental variable
  19.  *      (defined in the header file gmt_os2.h) to be used to specify the
  20.  *      data path for the GMT data files required.
  21.  *----------------------------------------------------------------------*/
  22.    char *gmt_data_path;
  23.    gmt_data_path = getenv (GMTDATAPATH);
  24.     sprintf (file, "%s/binned_shore_f.cdf\0", gmt_data_path);
  25.         fh[0] = !access (file, R_OK);
  26.     sprintf (file, "%s/binned_shore_h.cdf\0", gmt_data_path);
  27.         fh[1] = !access (file, R_OK);
  28.  
  29. /* -------------------Original coding below----------------------------
  30.  *
  31.  *    sprintf (file, "%s/binned_shore_f.cdf\0", LIBDIR);
  32.  *       fh[0] = !access (file, R_OK);
  33.  *
  34.  *    sprintf (file, "%s/binned_shore_h.cdf\0", LIBDIR);
  35.  *       fh[1] = !access (file, R_OK);
  36.  * --------------------------------------------------------------------*/
  37. }
  38.  
  39. int gmt_init_shore (res, c, w, e, s, n)
  40. char res;    /* Resolution (f, h, i, l, c */
  41. struct SHORE *c;
  42. double w, e, s, n; {
  43.     int i, nb, idiv, iw, ie, is, in, this_south, this_west;
  44.     short *stmp;
  45.     nclong *itmp;
  46.     nclong start[1], count[1];
  47.     char file[80];
  48. /*----------------------------------------------------------------------
  49.  *      Modifications made to permit the use of an environmental variable
  50.  *      (defined in the header file gmt_os2.h) to be used to specify the
  51.  *      data path for the GMT data files required.
  52.  *----------------------------------------------------------------------*/
  53.    char *gmt_data_path;
  54.    gmt_data_path = getenv (GMTDATAPATH);
  55.     sprintf (file, "%s/binned_shore_%c.cdf\0", gmt_data_path, res);
  56.     
  57. /* -------------------Original coding below----------------------------
  58.  *    sprintf (file, "%s/binned_shore_%c.cdf\0", LIBDIR, res);
  59.  * --------------------------------------------------------------------*/
  60.     
  61.         if (access (file, R_OK)) return (-1);
  62.         if ((c->cdfid = ncopen (file, NC_NOWRITE)) == -1) return (-1);
  63.                 
  64.     /* Get all id tags */
  65.     
  66.     c->bin_size_id    = ncvarid (c->cdfid, "Bin size in minutes");
  67.     c->bin_nx_id    = ncvarid (c->cdfid, "N bins in 0-360 longitude range");
  68.     c->bin_ny_id    = ncvarid (c->cdfid, "N bins in -90 +90 latitude range");
  69.     c->n_bin_id    = ncvarid (c->cdfid, "N bins in file");
  70.     c->n_seg_id    = ncvarid (c->cdfid, "N segments in file");
  71.     c->n_pt_id    = ncvarid (c->cdfid, "N points in file");
  72.     
  73.     c->bin_firstseg_id    = ncvarid (c->cdfid, "Id of first segment in a bin");
  74.     c->bin_info_id    = ncvarid (c->cdfid, "Embedded node levels in a bin");
  75.     c->bin_nseg_id    = ncvarid (c->cdfid, "N segments in a bin");
  76.     
  77.     c->seg_info_id    = ncvarid (c->cdfid, "Embedded <npts, levels, exit, entry> for a segment");
  78.     c->seg_area_id    = ncvarid (c->cdfid, "Area in 0.1km^2 of the parent polygon of a segment");
  79.     c->seg_start_id    = ncvarid (c->cdfid, "Id of first point in a segment");
  80.  
  81.     c->pt_dx_id    = ncvarid (c->cdfid, "Relative longitude from SE corner of bin");
  82.     c->pt_dy_id    = ncvarid (c->cdfid, "Relative latitude from SE corner of bin");
  83.  
  84.     /* Get attributes */
  85.  
  86.     ncattget (c->cdfid, c->pt_dx_id, "units", (void *)c->units);
  87.     ncattget (c->cdfid, NC_GLOBAL, "title", (void *)c->title);
  88.     ncattget (c->cdfid, NC_GLOBAL, "source", (void *)c->source);
  89.  
  90.     /* Get global variables */
  91.  
  92.     start[0] = 0;
  93.     ncvarget1 (c->cdfid, c->bin_size_id, start, (void *)&c->bin_size);
  94.     ncvarget1 (c->cdfid, c->bin_nx_id, start, (void *)&c->bin_nx);
  95.     ncvarget1 (c->cdfid, c->bin_ny_id, start, (void *)&c->bin_ny);
  96.     ncvarget1 (c->cdfid, c->n_bin_id, start, (void *)&c->n_bin);
  97.     ncvarget1 (c->cdfid, c->n_seg_id, start, (void *)&c->n_seg);
  98.     ncvarget1 (c->cdfid, c->n_pt_id, start, (void *)&c->n_pt);
  99.  
  100.     c->scale = (c->bin_size / 60.0) / 65535.0;
  101.     c->bsize = c->bin_size / 60.0;
  102.  
  103.     c->bins = (int *) memory ((char *)NULL, (int)c->n_bin, sizeof (int), "gmt_init_shore");
  104.     
  105.     /* Round off area to nearest multiple of block-dimension */
  106.         
  107.     iw = floor (w / c->bsize) * c->bsize;
  108.     ie =  ceil (e / c->bsize) * c->bsize;
  109.     is = 90 - ceil ((90.0 - s) / c->bsize) * c->bsize;
  110.     in = 90 - floor ((90.0 - n) / c->bsize) * c->bsize;
  111.     idiv = rint (360.0 / c->bsize);    /* Number of blocks per latitude band */
  112.     
  113.     for (i = nb = 0; i < c->n_bin; i++) {    /* Find which bins are needed */
  114.         this_south = 90 - c->bsize * ((i / idiv) + 1);
  115.         if (this_south < is || this_south >= in) continue;
  116.         this_west = c->bsize * (i % idiv) - 360;
  117.         while (this_west < iw) this_west += 360;
  118.         if (this_west >= ie) continue;
  119.         c->bins[nb] = i;
  120.         nb++;
  121.     }
  122.     c->bins = (int *) memory ((char *)c->bins, nb, sizeof (int), "gmt_init_shore");
  123.     c->nb = nb;
  124.     
  125.     /* Get bin variables, then extract only those corresponding to the bins to use */
  126.  
  127.     /* Allocate space for arrays of bin information */
  128.     
  129.     c->bin_info     = (short *) memory ((char *)NULL, nb, sizeof (short), "gmt_init_shore");
  130.     c->bin_nseg     = (short *) memory ((char *)NULL, nb, sizeof (short), "gmt_init_shore");
  131.     c->bin_firstseg     = (nclong *) memory ((char *)NULL, nb, sizeof (nclong), "gmt_init_shore");
  132.     
  133.     count[0] = c->n_bin;
  134.     stmp = (short *) memory ((char *)NULL, (int)c->n_bin, sizeof (short), "gmt_init_shore");
  135.     
  136.     ncvarget (c->cdfid, c->bin_info_id, start, count, (void *)stmp);
  137.     for (i = 0; i < c->nb; i++) c->bin_info[i] = stmp[c->bins[i]];
  138.     
  139.     ncvarget (c->cdfid, c->bin_nseg_id, start, count, (void *)stmp);
  140.     for (i = 0; i < c->nb; i++) c->bin_nseg[i] = stmp[c->bins[i]];
  141.     free ((char *)stmp);
  142.     
  143.     itmp = (nclong *) memory ((char *)NULL, (int)c->n_bin, sizeof (nclong), "gmt_init_shore");
  144.     ncvarget (c->cdfid, c->bin_firstseg_id, start, count, (void *)itmp);
  145.     for (i = 0; i < c->nb; i++) c->bin_firstseg[i] = itmp[c->bins[i]];
  146.     
  147.     free ((char *)itmp);
  148.     
  149.     return (0);
  150. }
  151.  
  152. int gmt_get_shore_bin (b, c, min_area, max_level)
  153. int b;            /* index number into c->bins */
  154. struct SHORE *c;
  155. double min_area;    /* Polygons with area less than this are ignored */
  156. int max_level; {    /* Polygons with higher levels are ignored */
  157.     nclong start[1], count[1];
  158.     nclong cut_area, *seg_area, *seg_info, *seg_start;
  159.     int s, i;
  160.     
  161.     c->node_level[0] = MIN (((ushort)c->bin_info[b] >> 9) & 7, max_level);
  162.     c->node_level[1] = MIN (((ushort)c->bin_info[b] >> 6) & 7, max_level);
  163.     c->node_level[2] = MIN (((ushort)c->bin_info[b] >> 3) & 7, max_level);
  164.     c->node_level[3] = MIN ((ushort)c->bin_info[b] & 7, max_level);
  165.     c->lon_sw = (c->bins[b] % c->bin_nx) * c->bin_size / 60.0;
  166.     c->lat_sw = 90.0 - ((c->bins[b] / c->bin_nx) + 1) * c->bin_size / 60.0;
  167.     c->ns = 0;
  168.     
  169.     if (c->bin_nseg[b] == 0) return;
  170.     
  171.     cut_area = rint (10.0 * min_area);
  172.     start[0] = c->bin_firstseg[b];
  173.     count[0] = c->bin_nseg[b];
  174.     
  175.     seg_area = (nclong *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (nclong), "gmt_get_shore_bin");
  176.     seg_info = (nclong *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (nclong), "gmt_get_shore_bin");
  177.     seg_start = (nclong *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (nclong), "gmt_get_shore_bin");
  178.     
  179.     ncvarget (c->cdfid, c->seg_area_id, start, count, (void *)seg_area);
  180.     ncvarget (c->cdfid, c->seg_info_id, start, count, (void *)seg_info);
  181.     ncvarget (c->cdfid, c->seg_start_id, start, count, (void *)seg_start);
  182.     
  183.     /* First tally how many useful segments */
  184.     
  185.     for (s = i = 0; i < c->bin_nseg[b]; i++) {
  186.         if (cut_area > 0 && seg_area[i] < cut_area) continue;
  187.         if (((seg_info[i] >> 6) & 7) > max_level) continue;
  188.         seg_area[s] = seg_area[i];
  189.         seg_info[s] = seg_info[i];
  190.         seg_start[s] = seg_start[i];
  191.         s++;
  192.     }
  193.     c->ns = s;
  194.     
  195.     if (c->ns == 0) {    /* No useful segments in this bin */
  196.         free ((char *) seg_info);    
  197.         free ((char *) seg_area);    
  198.         free ((char *) seg_start);
  199.         return;
  200.     }
  201.     
  202.     c->seg = (struct SHORE_SEGMENT *) memory ((char *)NULL, c->ns, sizeof (struct SHORE_SEGMENT), "gmt_get_shore_bin");
  203.     
  204.     for (s = 0; s < c->ns; s++) {
  205.         c->seg[s].level = (seg_info[s] >> 6) & 7;
  206.         c->seg[s].n = (seg_info[s] >> 9);
  207.         c->seg[s].entry = (seg_info[s] >> 3) & 7;
  208.         c->seg[s].exit = seg_info[s] & 7;
  209.         c->seg[s].dx = (short *) memory ((char *)NULL, (int)c->seg[s].n, sizeof (short), "gmt_get_shore_bin");
  210.         c->seg[s].dy = (short *) memory ((char *)NULL, (int)c->seg[s].n, sizeof (short), "gmt_get_shore_bin");
  211.         start[0] = seg_start[s];
  212.         count[0] = c->seg[s].n;
  213.         ncvarget (c->cdfid, c->pt_dx_id, start, count, (void *)c->seg[s].dx);
  214.         ncvarget (c->cdfid, c->pt_dy_id, start, count, (void *)c->seg[s].dy);
  215.             
  216.     }
  217.         
  218.     free ((char *) seg_info);    
  219.     free ((char *) seg_area);    
  220.     free ((char *) seg_start);    
  221.     
  222. }
  223.  
  224. int gmt_init_br (which, res, c, w, e, s, n)
  225. char which;    /* r(iver) or b(order) */
  226. char res;    /* Resolution (f, h, i, l, c */
  227. struct BR *c;
  228. double w, e, s, n; {
  229.     int i, nb, idiv, iw, ie, is, in, this_south, this_west;
  230.     short *stmp;
  231.     nclong *itmp;
  232.     nclong start[1], count[1];
  233.     char file[80];
  234.     char *gmt_data_path;
  235.    gmt_data_path = getenv (GMTDATAPATH);
  236.  
  237.     if (which == 'r')
  238. /*        sprintf (file, "%s/binned_river_%c.cdf\0", LIBDIR, res);   */
  239.         sprintf (file, "%s/binned_river_%c.cdf\0", gmt_data_path, res);
  240.     else
  241. /*        sprintf (file, "%s/binned_border_%c.cdf\0", LIBDIR, res);  */
  242.         sprintf (file, "%s/binned_border_%c.cdf\0", gmt_data_path, res);
  243.     
  244.         if (access (file, R_OK)) return (-1);
  245.         if ((c->cdfid = ncopen (file, NC_NOWRITE)) == -1) return (-1);
  246.         
  247.     /* Get all id tags */
  248.     
  249.     c->bin_size_id    = ncvarid (c->cdfid, "Bin size in minutes");
  250.     c->bin_nx_id    = ncvarid (c->cdfid, "N bins in 0-360 longitude range");
  251.     c->bin_ny_id    = ncvarid (c->cdfid, "N bins in -90 +90 latitude range");
  252.     c->n_bin_id    = ncvarid (c->cdfid, "N bins in file");
  253.     c->n_seg_id    = ncvarid (c->cdfid, "N segments in file");
  254.     c->n_pt_id    = ncvarid (c->cdfid, "N points in file");
  255.     
  256.     c->bin_firstseg_id    = ncvarid (c->cdfid, "Id of first segment in a bin");
  257.     c->bin_nseg_id    = ncvarid (c->cdfid, "N segments in a bin");
  258.     
  259.     c->seg_n_id    = ncvarid (c->cdfid, "N points for a segment");
  260.     c->seg_level_id    = ncvarid (c->cdfid, "Hierarchial level of a segment");
  261.     c->seg_start_id    = ncvarid (c->cdfid, "Id of first point in a segment");
  262.  
  263.     c->pt_dx_id    = ncvarid (c->cdfid, "Relative longitude from SE corner of bin");
  264.     c->pt_dy_id    = ncvarid (c->cdfid, "Relative latitude from SE corner of bin");
  265.  
  266.     /* Get attributes */
  267.  
  268.     ncattget (c->cdfid, c->pt_dx_id, "units", (void *)c->units);
  269.     ncattget (c->cdfid, NC_GLOBAL, "title", (void *)c->title);
  270.     ncattget (c->cdfid, NC_GLOBAL, "source", (void *)c->source);
  271.  
  272.     /* Get global variables */
  273.  
  274.     start[0] = 0;
  275.     ncvarget1 (c->cdfid, c->bin_size_id, start, (void *)&c->bin_size);
  276.     ncvarget1 (c->cdfid, c->bin_nx_id, start, (void *)&c->bin_nx);
  277.     ncvarget1 (c->cdfid, c->bin_ny_id, start, (void *)&c->bin_ny);
  278.     ncvarget1 (c->cdfid, c->n_bin_id, start, (void *)&c->n_bin);
  279.     ncvarget1 (c->cdfid, c->n_seg_id, start, (void *)&c->n_seg);
  280.     ncvarget1 (c->cdfid, c->n_pt_id, start, (void *)&c->n_pt);
  281.  
  282.     c->scale = (c->bin_size / 60.0) / 65535.0;
  283.     c->bsize = c->bin_size / 60.0;
  284.  
  285.     c->bins = (int *) memory ((char *)NULL, (int)c->n_bin, sizeof (int), "gmt_init_br");
  286.     
  287.     /* Round off area to nearest multiple of block-dimension */
  288.         
  289.     iw = floor (w / c->bsize) * c->bsize;
  290.     ie =  ceil (e / c->bsize) * c->bsize;
  291.     is = 90 - ceil ((90.0 - s) / c->bsize) * c->bsize;
  292.     in = 90 - floor ((90.0 - n) / c->bsize) * c->bsize;
  293.     idiv = rint (360.0 / c->bsize);    /* Number of blocks per latitude band */
  294.     
  295.     for (i = nb = 0; i < c->n_bin; i++) {    /* Find which bins are needed */
  296.         this_south = 90 - c->bsize * ((i / idiv) + 1);
  297.         if (this_south < is || this_south >= in) continue;
  298.         this_west = c->bsize * (i % idiv) - 360;
  299.         while (this_west < iw) this_west += 360;
  300.         if (this_west >= ie) continue;
  301.         c->bins[nb] = i;
  302.         nb++;
  303.     }
  304.     c->bins = (int *) memory ((char *)c->bins, nb, sizeof (int), "gmt_init_br");
  305.     c->nb = nb;
  306.     
  307.     /* Get bin variables, then extract only those corresponding to the bins to use */
  308.  
  309.     /* Allocate space for arrays of bin information */
  310.     
  311.     c->bin_nseg     = (short *) memory ((char *)NULL, nb, sizeof (short), "gmt_init_br");
  312.     c->bin_firstseg     = (nclong *) memory ((char *)NULL, nb, sizeof (nclong), "gmt_init_br");
  313.     
  314.     count[0] = c->n_bin;
  315.     stmp = (short *) memory ((char *)NULL, (int)c->n_bin, sizeof (short), "gmt_init_br");
  316.     
  317.     ncvarget (c->cdfid, c->bin_nseg_id, start, count, (void *)stmp);
  318.     for (i = 0; i < c->nb; i++) c->bin_nseg[i] = stmp[c->bins[i]];
  319.     free ((char *)stmp);
  320.     
  321.     itmp = (nclong *) memory ((char *)NULL, (int)c->n_bin, sizeof (nclong), "gmt_init_br");
  322.     ncvarget (c->cdfid, c->bin_firstseg_id, start, count, (void *)itmp);
  323.     for (i = 0; i < c->nb; i++) c->bin_firstseg[i] = itmp[c->bins[i]];
  324.     
  325.     free ((char *)itmp);
  326.     
  327.     return (0);
  328. }
  329.  
  330. int gmt_get_br_bin (b, c, level, n_levels)
  331. int b;            /* index number into c->bins */
  332. struct BR *c;
  333. int level[];        /* Levels of features to extract */
  334. int n_levels;        /* # of such levels. 0 means use all levels */
  335. {
  336.     nclong start[1], count[1];
  337.     nclong *seg_start;
  338.     short *seg_n, *seg_level;
  339.     int s, i, k, skip;
  340.     
  341.     c->lon_sw = (c->bins[b] % c->bin_nx) * c->bin_size / 60.0;
  342.     c->lat_sw = 90.0 - ((c->bins[b] / c->bin_nx) + 1) * c->bin_size / 60.0;
  343.     c->ns = c->bin_nseg[b];
  344.     
  345.     if (c->ns == 0) return;
  346.     
  347.     start[0] = c->bin_firstseg[b];
  348.     count[0] = c->bin_nseg[b];
  349.     
  350.     seg_n = (short *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (short), "gmt_get_br_bin");
  351.     seg_level = (short *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (short), "gmt_get_br_bin");
  352.     seg_start = (nclong *) memory ((char *)NULL, (int)c->bin_nseg[b], sizeof (nclong), "gmt_get_br_bin");
  353.     
  354.     ncvarget (c->cdfid, c->seg_n_id, start, count, (void *)seg_n);
  355.     ncvarget (c->cdfid, c->seg_level_id, start, count, (void *)seg_level);
  356.     ncvarget (c->cdfid, c->seg_start_id, start, count, (void *)seg_start);
  357.     
  358.     c->seg = (struct BR_SEGMENT *) memory ((char *)NULL, c->ns, sizeof (struct BR_SEGMENT), "gmt_get_br_bin");
  359.     
  360.     for (s = i = 0; i < c->ns; i++) {
  361.         if (n_levels == 0)
  362.             skip = FALSE;
  363.         else {
  364.             for (k = 0, skip = TRUE; skip && k < n_levels; k++)
  365.                 if (seg_level[i] == level[k]) skip = FALSE;
  366.         }
  367.         if (skip) continue;
  368.         
  369.         c->seg[s].n = seg_n[i];
  370.         c->seg[s].level = seg_level[i];
  371.         c->seg[s].dx = (short *) memory ((char *)NULL, (int)c->seg[s].n, sizeof (short), "gmt_get_br_bin");
  372.         c->seg[s].dy = (short *) memory ((char *)NULL, (int)c->seg[s].n, sizeof (short), "gmt_get_br_bin");
  373.         start[0] = seg_start[i];
  374.         count[0] = c->seg[s].n;
  375.         ncvarget (c->cdfid, c->pt_dx_id, start, count, (void *)c->seg[s].dx);
  376.         ncvarget (c->cdfid, c->pt_dy_id, start, count, (void *)c->seg[s].dy);
  377.         
  378.         s++;
  379.     }
  380.     
  381.     c->ns = s;
  382.  
  383.     free ((char *) seg_n);    
  384.     free ((char *) seg_level);    
  385.     free ((char *) seg_start);    
  386.     
  387. }
  388.  
  389. int get_position (side, x, y)
  390. int side;
  391. short x, y; {
  392.     /* Returns the position along the given side */
  393.     
  394.     return ((side%2) ? ((side == 1) ? (ushort)y : MAX_DELTA - (ushort)y) : ((side == 0) ? (ushort)x : MAX_DELTA - (ushort)x));
  395. }
  396.  
  397. int gmt_prepare_sides (c, dir)
  398. struct SHORE *c;
  399. int dir;
  400. {
  401.     int s, i, n[4], asc_sort(), desc_sort();
  402.     
  403.     c->lon_corner[0] = c->lon_sw + ((dir == 1) ? c->bsize : 0.0);
  404.     c->lon_corner[1] = c->lon_sw + c->bsize;
  405.     c->lon_corner[2] = c->lon_sw + ((dir == 1) ? 0.0 : c->bsize);
  406.     c->lon_corner[3] = c->lon_sw;
  407.     c->lat_corner[0] = c->lat_sw;
  408.     c->lat_corner[1] = c->lat_sw + ((dir == 1) ? c->bsize : 0.0);
  409.     c->lat_corner[2] = c->lat_sw + c->bsize;
  410.     c->lat_corner[3] = c->lat_sw + ((dir == 1) ? 0.0 : c->bsize);
  411.  
  412.     for (i = 0; i < 4; i++) c->nside[i] = n[i] = 1;
  413.     /* for (s = 0; s < c->ns; s++) if (c->seg[s].level < 3 && c->seg[s].entry < 4) c->nside[c->seg[s].entry]++; */
  414.     for (s = 0; s < c->ns; s++) if (c->seg[s].entry < 4) c->nside[c->seg[s].entry]++;
  415.     
  416.     for (i = c->n_entries = 0; i < 4; i++) {    /* Allocate memory and add corners */
  417.         c->side[i] = (struct SIDE *) memory (CNULL, c->nside[i], sizeof (struct SIDE), "gmt_prepare_sides");
  418.         c->side[i][0].pos = (dir == 1) ? MAX_DELTA : 0;
  419.         c->side[i][0].id = i - 4;
  420.         c->n_entries += c->nside[i] - 1;
  421.     }
  422.         
  423.     for (s = 0; s < c->ns; s++) {    /* Add entry points */
  424.         /* if (c->seg[s].level > 2 || (i = c->seg[s].entry) == 4) continue; */
  425.         if ((i = c->seg[s].entry) == 4) continue;
  426.         c->side[i][n[i]].pos = get_position (i, c->seg[s].dx[0], c->seg[s].dy[0]);
  427.         c->side[i][n[i]].id = s;
  428.         n[i]++;
  429.     }
  430.     
  431.     for (i = 0; i < 4; i++)    {    /* sort on position */
  432.         if (dir == 1)
  433.             qsort ((char *)c->side[i], c->nside[i], sizeof (struct SIDE), asc_sort);
  434.         else
  435.             qsort ((char *)c->side[i], c->nside[i], sizeof (struct SIDE), desc_sort);
  436.     }
  437. }
  438.  
  439. int gmt_assemble_shore (c, dir, assemble, shift, edge, pol)
  440. struct SHORE *c;
  441. int dir;
  442. BOOLEAN assemble;    /* TRUE if polygons is needed */
  443. BOOLEAN shift;        /* TRUE if longitudes may have to be shifted */
  444. double edge;        /* Edge test for shifting */
  445. struct POL *pol[];
  446. {
  447.     struct POL *p;
  448.     int start_side, next_side, id, P = 0, more, p_alloc, i0, i1, j, wet_or_dry, use_this_level;
  449.     int n_alloc, cid, nid, add, first_pos, entry_pos, n, low_level, high_level;
  450.     BOOLEAN completely_inside;
  451.     double *xtmp, *ytmp, plon, plat;
  452.     
  453.     if (!assemble) {    /* Easy, just need to scale all segments to degrees and return */
  454.     
  455.         p = (struct POL *) memory (CNULL, c->ns, sizeof (struct POL), "gmt_assemble_shore");
  456.         
  457.         for (id = 0; id < c->ns; id++) {
  458.             p[id].lon = (double *) memory (CNULL, (int)c->seg[id].n, sizeof (double), "gmt_assemble_shore");
  459.             p[id].lat = (double *) memory (CNULL, (int)c->seg[id].n, sizeof (double), "gmt_assemble_shore");
  460.             p[id].n = copy_to_shore_path (p[id].lon, p[id].lat, c, id);
  461.             p[id].level = c->seg[id].level;
  462.             p[id].interior = FALSE;
  463.             if (shift) gmt_path_shift (p[id].lon, p[id].lat, p[id].n, edge);
  464.         }
  465.     
  466.         *pol = p;
  467.         return (c->ns);
  468.     }
  469.     
  470.     for (n = high_level = 0; n < 4; n++) high_level = MAX (c->node_level[n], high_level);
  471.     wet_or_dry = (dir == 1) ? 1 : 0;
  472.     use_this_level = (high_level%2 == wet_or_dry);
  473.  
  474.     if (c->ns == 0 && !use_this_level) return (0);    /* No polygons for this bin */
  475.     
  476.     /* Here we must assemble [at least one] polygons in the correct order */
  477.     
  478.     for (n = 0, completely_inside = TRUE; completely_inside && n < c->ns; n++) if (c->seg[n].entry != 4) completely_inside = FALSE;
  479.  
  480.     gmt_prepare_sides (c, dir);
  481.     
  482.     p_alloc = (c->ns == 0) ? 1 : GMT_SMALL_CHUNK;
  483.     p = (struct POL *) memory (CNULL, p_alloc, sizeof (struct POL), "gmt_assemble_shore");
  484.     
  485.     low_level = MAX_LEVEL;
  486.     
  487.     if (completely_inside && use_this_level) {    /* Must include path of this bin outline as first polygon */
  488.         for (i0 = j = n = n_alloc = 0; j < 4; j++) {
  489.             i1 = (i0 + dir + 4) % 4;
  490.             if ((add = map_path (c->lon_corner[i0], c->lat_corner[i0], c->lon_corner[i1], c->lat_corner[i1], &xtmp, &ytmp))) {
  491.                 n_alloc += add;
  492.                 p[0].lon = (double *) memory ((char *)p[0].lon, n_alloc, sizeof (double), "gmt_assemble_shore");
  493.                 p[0].lat = (double *) memory ((char *)p[0].lat, n_alloc, sizeof (double), "gmt_assemble_shore");
  494.                 memcpy ((char *)&p[0].lon[n], (char *)xtmp, add * sizeof (double));
  495.                 memcpy ((char *)&p[0].lat[n], (char *)ytmp, add * sizeof (double));
  496.                 n += add;
  497.                 free ((char *)xtmp);
  498.                 free ((char *)ytmp);
  499.             }
  500.             i0 = i1;
  501.         }
  502.         p[0].n = n;
  503.         p[0].level = c->node_level[0];    /* Any corner will do */
  504.         p[0].interior = FALSE;
  505.         P = 1;
  506.     }        
  507.     
  508.     while (c->n_entries > 0) {    /* More segments to connect */
  509.     
  510.         low_level = MAX_LEVEL;
  511.         start_side = 0;
  512.         id = get_first_entry (c, dir, &start_side);
  513.         next_side = c->seg[id].exit;
  514.         
  515.         n_alloc = c->seg[id].n;
  516.         p[P].lon = (double *) memory (CNULL, n_alloc, sizeof (double), "gmt_assemble_shore");
  517.         p[P].lat = (double *) memory (CNULL, n_alloc, sizeof (double), "gmt_assemble_shore");
  518.         n = copy_to_shore_path (p[P].lon, p[P].lat, c, id);
  519.         if (c->seg[id].level < low_level) low_level = c->seg[id].level;
  520.         
  521.         more = TRUE;
  522.         first_pos = get_position (start_side, c->seg[id].dx[0], c->seg[id].dy[0]);
  523.         while (more) {
  524.         
  525.             id = get_next_entry (c, dir, next_side, id);
  526.         
  527.             if (id < 0) {    /* Corner */
  528.                 cid = id + 4;
  529.                 nid = (dir == 1) ? (cid + 1) % 4 : cid;
  530.                 if ((add = map_path (p[P].lon[n-1], p[P].lat[n-1], c->lon_corner[cid], c->lat_corner[cid], &xtmp, &ytmp))) {
  531.                     n_alloc += add;
  532.                     p[P].lon = (double *) memory ((char *)p[P].lon, n_alloc, sizeof (double), "gmt_assemble_shore");
  533.                     p[P].lat = (double *) memory ((char *)p[P].lat, n_alloc, sizeof (double), "gmt_assemble_shore");
  534.                     memcpy ((char *)&p[P].lon[n], (char *)xtmp, add * sizeof (double));
  535.                     memcpy ((char *)&p[P].lat[n], (char *)ytmp, add * sizeof (double));
  536.                     n += add;
  537.                 }
  538.                 next_side = ((id + 4) + dir + 4) % 4;
  539.                 if (c->node_level[nid] < low_level) low_level = c->node_level[nid];
  540.             }
  541.             else {
  542.                 shore_to_degree (c, c->seg[id].dx[0], c->seg[id].dy[0], &plon, &plat);
  543.                 if ((add = map_path (p[P].lon[n-1], p[P].lat[n-1], plon, plat, &xtmp, &ytmp))) {
  544.                     n_alloc += add;
  545.                     p[P].lon = (double *) memory ((char *)p[P].lon, n_alloc, sizeof (double), "gmt_assemble_shore");
  546.                     p[P].lat = (double *) memory ((char *)p[P].lat, n_alloc, sizeof (double), "gmt_assemble_shore");
  547.                     memcpy ((char *)&p[P].lon[n], (char *)xtmp, add * sizeof (double));
  548.                     memcpy ((char *)&p[P].lat[n], (char *)ytmp, add * sizeof (double));
  549.                     n += add;
  550.                 }
  551.                 entry_pos = get_position (next_side, c->seg[id].dx[0], c->seg[id].dy[0]);
  552.                 if (next_side == start_side && entry_pos == first_pos)
  553.                     more = FALSE;
  554.                 else {
  555.                     n_alloc += c->seg[id].n;
  556.                     p[P].lon = (double *) memory ((char *)p[P].lon, n_alloc, sizeof (double), "gmt_assemble_shore");
  557.                     p[P].lat = (double *) memory ((char *)p[P].lat, n_alloc, sizeof (double), "gmt_assemble_shore");
  558.                     n += copy_to_shore_path (&p[P].lon[n], &p[P].lat[n], c, id);
  559.                     next_side = c->seg[id].exit;
  560.                     if (c->seg[id].level < low_level) low_level = c->seg[id].level;
  561.                 }
  562.             }
  563.             if (add) {
  564.                 free ((char *)xtmp);
  565.                 free ((char *)ytmp);
  566.             }
  567.         }
  568.         p[P].n = n;
  569.         p[P].interior = FALSE;
  570.         p[P].level = (dir == 1) ? 2 * ((low_level - 1) / 2) + 1: 2 * (low_level/2);
  571.         P++;
  572.         if (P == p_alloc) {
  573.             p_alloc += GMT_SMALL_CHUNK;
  574.             p = (struct POL *) memory ((char *)p, p_alloc, sizeof (struct POL), "gmt_assemble_shore");
  575.         }
  576.         
  577.     }
  578.     
  579.     /* Then add all interior polygons, if any */
  580.     
  581.     for (id = 0; id < c->ns; id++) {
  582.         if (c->seg[id].entry < 4) continue;
  583.         n_alloc = c->seg[id].n;
  584.         p[P].lon = (double *) memory (CNULL, n_alloc, sizeof (double), "gmt_assemble_shore");
  585.         p[P].lat = (double *) memory (CNULL, n_alloc, sizeof (double), "gmt_assemble_shore");
  586.         p[P].n = copy_to_shore_path (p[P].lon, p[P].lat, c, id);
  587.         p[P].interior = TRUE;
  588.         p[P].level = c->seg[id].level;
  589.         P++;
  590.         if (P == p_alloc) {
  591.             p_alloc += GMT_SMALL_CHUNK;
  592.             p = (struct POL *) memory ((char *)p, p_alloc, sizeof (struct POL), "gmt_assemble_shore");
  593.         }
  594.     }
  595.     
  596.     gmt_pau_sides (c);
  597.  
  598.     if (c->ns > 0) p = (struct POL *) memory ((char *)p, P, sizeof (struct POL), "gmt_assemble_shore");
  599.     
  600.     if (shift) for (id = 0; id < P; id++) gmt_path_shift (p[id].lon, p[id].lat, p[id].n, edge);
  601.  
  602.     *pol = p;
  603.     return (P);
  604. }
  605.         
  606. int gmt_assemble_br (c, shift, edge, pol)
  607. struct BR *c;
  608. BOOLEAN shift;        /* TRUE if longitudes may have to be shifted */
  609. double edge;        /* Edge test for shifting */
  610. struct POL *pol[];
  611. {
  612.     struct POL *p;
  613.     int id;
  614.     
  615.     p = (struct POL *) memory (CNULL, c->ns, sizeof (struct POL), "gmt_assemble_br");
  616.     
  617.     for (id = 0; id < c->ns; id++) {
  618.         p[id].lon = (double *) memory (CNULL, (int)c->seg[id].n, sizeof (double), "gmt_assemble_br");
  619.         p[id].lat = (double *) memory (CNULL, (int)c->seg[id].n, sizeof (double), "gmt_assemble_br");
  620.         p[id].n = copy_to_br_path (p[id].lon, p[id].lat, c, id);
  621.         p[id].level = c->seg[id].level;
  622.         if (shift) gmt_path_shift (p[id].lon, p[id].lat, p[id].n, edge);
  623.     }
  624.     
  625.     *pol = p;
  626.     return (c->ns);
  627. }
  628.         
  629. int copy_to_shore_path (lon, lat, s, id)
  630. double lon[];
  631. double lat[];
  632. struct SHORE *s;
  633. int id; {
  634.     int i;
  635.     for (i = 0; i < s->seg[id].n; i++)
  636.         shore_to_degree (s, s->seg[id].dx[i], s->seg[id].dy[i], &lon[i], &lat[i]);
  637.     return (s->seg[id].n);
  638. }
  639.  
  640. int copy_to_br_path (lon, lat, s, id)
  641. double lon[];
  642. double lat[];
  643. struct BR *s;
  644. int id; {
  645.     int i;
  646.     for (i = 0; i < s->seg[id].n; i++)
  647.         br_to_degree (s, s->seg[id].dx[i], s->seg[id].dy[i], &lon[i], &lat[i]);
  648.     return (s->seg[id].n);
  649. }
  650.  
  651. int shore_to_degree (c, dx, dy, lon, lat)
  652. struct SHORE *c;
  653. short dx, dy;
  654. double *lon, *lat; {
  655.     *lon = c->lon_sw + ((ushort)dx) * c->scale;
  656.     *lat = c->lat_sw + ((ushort)dy) * c->scale;
  657. }
  658.  
  659. int br_to_degree (c, dx, dy, lon, lat)
  660. struct BR *c;
  661. short dx, dy;
  662. double *lon, *lat; {
  663.     *lon = c->lon_sw + ((ushort)dx) * c->scale;
  664.     *lat = c->lat_sw + ((ushort)dy) * c->scale;
  665. }
  666.  
  667. int get_next_entry (c, dir, side, id)
  668. struct SHORE *c;
  669. int dir, id;
  670. int side; {
  671.     /* Finds the next entry point on the given side that is further away
  672.      * in the <dir> direction than previous point.  It removes the info
  673.      * regarding the new entry from the SIDE structure */
  674.      
  675.     int k, pos, n;
  676.     
  677.     if (id < 0)
  678.         pos = (dir == 1) ? 0 : MAX_DELTA;
  679.     else {
  680.         n = c->seg[id].n - 1;
  681.         pos = get_position (side, c->seg[id].dx[n], c->seg[id].dy[n]);
  682.     }
  683.  
  684.     if (dir == 1) {
  685.         for (k = 0; k < c->nside[side] && c->side[side][k].pos < pos; k++);
  686.         id = c->side[side][k].id;
  687.         for (k++; k < c->nside[side]; k++) c->side[side][k-1] = c->side[side][k];
  688.         c->nside[side]--;
  689.     }
  690.     else {
  691.         for (k = 0; k < c->nside[side] && c->side[side][k].pos > pos; k++);
  692.         id = c->side[side][k].id;
  693.         for (k++; k < c->nside[side]; k++) c->side[side][k-1] = c->side[side][k];
  694.         c->nside[side]--;
  695.     }
  696.     if (id >= 0) c->n_entries--;
  697.     return (id);
  698. }
  699.  
  700. int get_first_entry (c, dir, side)
  701. struct SHORE *c;
  702. int dir;
  703. int *side; {
  704.     int try = 0;
  705.     while (try < 4 && (c->nside[*side] == 0 || (c->nside[*side] == 1 && c->side[*side][0].id < 0))) {    /* No entries or only a corner left on this side */
  706.         try++;
  707.         *side = (*side + dir + 4) % 4;
  708.     }
  709.     if (try == 4) return (-5);
  710.     return (c->side[*side][0].id);
  711. }    
  712.  
  713. int asc_sort (a, b)
  714. struct SIDE *a, *b; {
  715.     if (a->pos < b->pos) return (-1);
  716.     if (a->pos > b->pos) return (1);
  717.     return (0);
  718. }
  719.  
  720. int desc_sort (a, b)
  721. struct SIDE *a, *b; {
  722.     if (a->pos < b->pos) return (1);
  723.     if (a->pos > b->pos) return (-1);
  724.     return (0);
  725. }
  726.  
  727. int free_shore (c)
  728. struct SHORE *c; {    /* Removes allocated variables for this block only */
  729.     int i;
  730.     
  731.     for (i = 0; i < c->ns; i++) {
  732.         free ((char *)c->seg[i].dx);
  733.         free ((char *)c->seg[i].dy);
  734.     }
  735.     if (c->ns) free ((char *)c->seg);
  736. }
  737.         
  738. int free_br (c)
  739. struct BR *c; {    /* Removes allocated variables for this block only */
  740.     int i;
  741.     
  742.     for (i = 0; i < c->ns; i++) {
  743.         free ((char *)c->seg[i].dx);
  744.         free ((char *)c->seg[i].dy);
  745.     }
  746.     if (c->ns) free ((char *)c->seg);
  747. }
  748.         
  749. int shore_cleanup (c)
  750. struct SHORE *c; {
  751.     free ((char *)c->bins);
  752.     free ((char *)c->bin_info);
  753.     free ((char *)c->bin_nseg);
  754.     free ((char *)c->bin_firstseg);
  755.     ncclose (c->cdfid);
  756. }
  757.  
  758. int br_cleanup (c)
  759. struct BR *c; {
  760.     free ((char *)c->bins);
  761.     free ((char *)c->bin_nseg);
  762.     free ((char *)c->bin_firstseg);
  763.     ncclose (c->cdfid);
  764. }
  765.  
  766. int gmt_pau_sides (c)
  767. struct SHORE *c;
  768. {
  769.     int i;
  770.     for (i = 0; i < 4; i++) free ((char *)c->side[i]);
  771. }
  772.  
  773. int free_polygons (p, n)
  774. struct POL p[];
  775. int n; {
  776.     int k;
  777.     for (k = 0; k < n; k++) {
  778.         free ((char *)p[k].lon);
  779.         free ((char *)p[k].lat);
  780.     }
  781. }
  782.  
  783. int gmt_path_shift (lon, lat, n, edge)
  784. double lon[], lat[];
  785. int n;
  786. double edge; {
  787.     int i;
  788.     
  789.     for (i = 0; i < n; i++) if (lon[i] >= edge) lon[i] -= 360.0;
  790. }
  791.  
  792.