/*--------------------------------------------------------------------
* The GMT-system: @(#)gmt_shore.c 3.78 12/01/99
*
* Copyright (c) 1991-1999 by P. Wessel and W. H. F. Smith
* See COPYING file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; version 2 of the License.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* Contact info: www.soest.hawaii.edu/gmt
*--------------------------------------------------------------------*/
#include "gmt.h"
/*
* These functions simplifies the access to the GMT shoreline, border, and river
* databases.
*
* The PUBLIC functions are:
*
* GMT_init_shore : Opens selected shoreline database and initializes structures
* GMT_get_shore_bin : Returns all selected shore data for this bin
* GMT_init_br : Opens selected border/river database and initializes structures
* GMT_get_br_bin : Returns all selected border/river data for this bin
* GMT_assemble_shore : Creates polygons or lines from shoreline segments
* GMT_prep_polygons : Wraps polygons if necessary and prepares them for use
* GMT_assemble_br : Creates lines from border or river segments
* GMT_free_shore : Frees up memory used by shorelines for this bin
* GMT_free_br : Frees up memory used by shorelines for this bin
* GMT_shore_cleanup : Frees up main shoreline structure memory
* GMT_br_cleanup : Frees up main river/border structure memory
*
* Author: Paul Wessel
* Date: 13-JUN-1995
* Modified: 16-AUG-1999
*
*/
int GMT_copy_to_shore_path (double *lon, double *lat, struct GMT_SHORE *s, int id);
int GMT_shore_get_first_entry (struct GMT_SHORE *c, int dir, int *side);
int GMT_shore_get_position (int side, short int x, short int y);
int GMT_shore_get_next_entry (struct GMT_SHORE *c, int dir, int side, int id);
int GMT_copy_to_br_path (double *lon, double *lat, struct GMT_BR *s, int id);
void GMT_shore_to_degree (struct GMT_SHORE *c, short int dx, short int dy, double *lon, double *lat);
void GMT_shore_pau_sides (struct GMT_SHORE *c);
void GMT_shore_path_shift (double *lon, double *lat, int n, double edge);
void GMT_shore_path_shift2 (double *lon, double *lat, int n, double west, double east, BOOLEAN leftmost);
void GMT_br_to_degree (struct GMT_BR *c, short int dx, short int dy, double *lon, double *lat);
void shore_prepare_sides(struct GMT_SHORE *c, int dir);
int GMT_shore_asc_sort (const void *a, const void *b);
int GMT_shore_desc_sort(const void *a, const void *b);
int check_nc_status (int status)
{
if (status != NC_NOERR) {
fprintf (stderr, "gmt_nan_init: %s\n", nc_strerror(status));
return (-1);
}
return 0;
}
int GMT_set_resolution (char *res, char opt)
{
/* Decodes the -D<res> option and returns the base integer value */
int base;
switch (*res) {
case 'f': /* Full */
base = 0;
break;
case 'h': /* High */
base = 1;
break;
case 'i': /* Intermediate */
base = 2;
break;
case 'l': /* Low */
base = 3;
break;
case 'c': /* Crude */
base = 4;
break;
default:
fprintf (stderr, "%s: GMT SYNTAX ERROR -%c option: Unknown modifier %c [Defaults to -%cl]\n", GMT_program, opt, *res, opt);
base = 3;
*res = 'l';
break;
}
return (base);
}
int GMT_init_shore (char res, struct GMT_SHORE *c, double w, double e, double s, double n)
/* res: Resolution (f, h, i, l, c */
{
int i, nb, idiv, iw, ie, is, in, this_south, this_west;
short *stmp;
int *itmp;
size_t start[1], count[1];
char file[32], path[BUFSIZ];
sprintf (file, "binned_GSHHS_%c.cdf\0", res);
if (!GMT_getpathname (file, path)) return (-1); /* Failed to find file */
check_nc_status (nc_open (path, NC_NOWRITE,&c->cdfid));
/* Get all id tags */
check_nc_status (nc_inq_varid (c->cdfid, "Bin_size_in_minutes", &c->bin_size_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_360_longitude_range", &c->bin_nx_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_180_degree_latitude_range", &c->bin_ny_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_file", &c->n_bin_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_segments_in_file", &c->n_seg_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_points_in_file", &c->n_pt_id));
check_nc_status (nc_inq_varid (c->cdfid, "Id_of_first_segment_in_a_bin", &c->bin_firstseg_id));
check_nc_status (nc_inq_varid (c->cdfid, "Embedded_node_levels_in_a_bin", &c->bin_info_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_segments_in_a_bin", &c->bin_nseg_id));
check_nc_status (nc_inq_varid (c->cdfid, "Embedded_npts_levels_exit_entry_for_a_segment", &c->seg_info_id));
check_nc_status (nc_inq_varid (c->cdfid, "Ten_times_the_km_squared_area_of_the_parent_polygon_of_a_segment", &c->seg_area_id));
check_nc_status (nc_inq_varid (c->cdfid, "Id_of_first_point_in_a_segment", &c->seg_start_id));
check_nc_status (nc_inq_varid (c->cdfid, "Relative_longitude_from_SW_corner_of_bin", &c->pt_dx_id));
check_nc_status (nc_inq_varid (c->cdfid, "Relative_latitude_from_SW_corner_of_bin", &c->pt_dy_id));
/* Get attributes */
check_nc_status (nc_get_att_text (c->cdfid, c->pt_dx_id, "units", c->units));
check_nc_status (nc_get_att_text (c->cdfid, NC_GLOBAL, "title", c->title));
check_nc_status (nc_get_att_text (c->cdfid, NC_GLOBAL, "source", c->source));
/* Get global variables */
start[0] = 0;
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_size_id, start, &c->bin_size));
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_nx_id, start, &c->bin_nx));
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_ny_id, start, &c->bin_ny));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_bin_id, start, &c->n_bin));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_seg_id, start, &c->n_seg));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_pt_id, start, &c->n_pt));
c->scale = (c->bin_size / 60.0) / 65535.0;
c->bsize = c->bin_size / 60.0;
c->bins = (int *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (int), "GMT_init_shore");
/* Round off area to nearest multiple of block-dimension */
iw = (int)(floor (w / c->bsize) * c->bsize);
ie = (int)(ceil (e / c->bsize) * c->bsize);
is = 90 - (int)(ceil ((90.0 - s) / c->bsize) * c->bsize);
in = 90 - (int)(floor ((90.0 - n) / c->bsize) * c->bsize);
idiv = irint (360.0 / c->bsize); /* Number of blocks per latitude band */
for (i = nb = 0; i < c->n_bin; i++) { /* Find which bins are needed */
this_south = 90 - (int)(c->bsize * ((i / idiv) + 1));
if (this_south < is || this_south >= in) continue;
this_west = (int)(c->bsize * (i % idiv)) - 360;
while (this_west < iw) this_west += 360;
if (this_west >= ie) continue;
c->bins[nb] = i;
nb++;
}
c->bins = (int *) GMT_memory ((void *)c->bins, (size_t)nb, sizeof (int), "GMT_init_shore");
c->nb = nb;
/* Get bin variables, then extract only those corresponding to the bins to use */
/* Allocate space for arrays of bin information */
c->bin_info = (short *) GMT_memory (VNULL, (size_t)nb, sizeof (short), "GMT_init_shore");
c->bin_nseg = (short *) GMT_memory (VNULL, (size_t)nb, sizeof (short), "GMT_init_shore");
c->bin_firstseg = (int *) GMT_memory (VNULL, (size_t)nb, sizeof (int), "GMT_init_shore");
count[0] = c->n_bin;
stmp = (short *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (short), "GMT_init_shore");
check_nc_status (nc_get_vara_short (c->cdfid, c->bin_info_id, start, count, stmp));
for (i = 0; i < c->nb; i++) c->bin_info[i] = stmp[c->bins[i]];
check_nc_status (nc_get_vara_short (c->cdfid, c->bin_nseg_id, start, count, stmp));
for (i = 0; i < c->nb; i++) c->bin_nseg[i] = stmp[c->bins[i]];
GMT_free ((void *)stmp);
itmp = (int *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (int), "GMT_init_shore");
check_nc_status (nc_get_vara_int (c->cdfid, c->bin_firstseg_id, start, count, itmp));
for (i = 0; i < c->nb; i++) c->bin_firstseg[i] = itmp[c->bins[i]];
GMT_free ((void *)itmp);
return (0);
}
void GMT_get_shore_bin (int b, struct GMT_SHORE *c, double min_area, int min_level, int max_level)
/* b: index number into c->bins */
/* min_area: Polygons with area less than this are ignored */
/* min_level: Polygons with lower levels are ignored */
/* max_level: Polygons with higher levels are ignored */
{
size_t start[1], count[1];
int cut_area, *seg_area, *seg_info, *seg_start;
int s, i;
double w, e, dx;
c->node_level[0] = (unsigned char)MIN (((unsigned short)c->bin_info[b] >> 9) & 7, max_level);
c->node_level[1] = (unsigned char)MIN (((unsigned short)c->bin_info[b] >> 6) & 7, max_level);
c->node_level[2] = (unsigned char)MIN (((unsigned short)c->bin_info[b] >> 3) & 7, max_level);
c->node_level[3] = (unsigned char)MIN ((unsigned short)c->bin_info[b] & 7, max_level);
dx = c->bin_size / 60.0;
c->lon_sw = (c->bins[b] % c->bin_nx) * dx;
c->lat_sw = 90.0 - ((c->bins[b] / c->bin_nx) + 1) * dx;
c->ns = 0;
/* Determine if this bin is one of the bins at the left side of the map */
w = c->lon_sw;
while (w > project_info.w && GMT_world_map) w -= 360.0;
e = w + dx;
c->leftmost_bin = ((w <= project_info.w) && (e > project_info.w));
if (c->bin_nseg[b] == 0) return;
cut_area = irint (10.0 * min_area);
start[0] = c->bin_firstseg[b];
count[0] = c->bin_nseg[b];
seg_area = (int *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (int), "GMT_get_shore_bin");
seg_info = (int *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (int), "GMT_get_shore_bin");
seg_start = (int *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (int), "GMT_get_shore_bin");
check_nc_status (nc_get_vara_int (c->cdfid, c->seg_area_id, start, count, seg_area));
check_nc_status (nc_get_vara_int (c->cdfid, c->seg_info_id, start, count, seg_info));
check_nc_status (nc_get_vara_int (c->cdfid, c->seg_start_id, start, count, seg_start));
/* First tally how many useful segments */
for (s = i = 0; i < c->bin_nseg[b]; i++) {
if (cut_area > 0 && seg_area[i] < cut_area) continue;
if (((seg_info[i] >> 6) & 7) < min_level) continue;
if (((seg_info[i] >> 6) & 7) > max_level) continue;
seg_area[s] = seg_area[i];
seg_info[s] = seg_info[i];
seg_start[s] = seg_start[i];
s++;
}
c->ns = s;
if (c->ns == 0) { /* No useful segments in this bin */
GMT_free ((void *) seg_info);
GMT_free ((void *) seg_area);
GMT_free ((void *) seg_start);
return;
}
c->seg = (struct GMT_SHORE_SEGMENT *) GMT_memory (VNULL, (size_t)c->ns, sizeof (struct GMT_SHORE_SEGMENT), "GMT_get_shore_bin");
for (s = 0; s < c->ns; s++) {
c->seg[s].level = (seg_info[s] >> 6) & 7;
c->seg[s].n = (seg_info[s] >> 9);
c->seg[s].entry = (seg_info[s] >> 3) & 7;
c->seg[s].exit = seg_info[s] & 7;
c->seg[s].dx = (short *) GMT_memory (VNULL, (size_t)c->seg[s].n, sizeof (short), "GMT_get_shore_bin");
c->seg[s].dy = (short *) GMT_memory (VNULL, (size_t)c->seg[s].n, sizeof (short), "GMT_get_shore_bin");
start[0] = seg_start[s];
count[0] = c->seg[s].n;
check_nc_status (nc_get_vara_short (c->cdfid, c->pt_dx_id, start, count, c->seg[s].dx));
check_nc_status (nc_get_vara_short (c->cdfid, c->pt_dy_id, start, count, c->seg[s].dy));
}
GMT_free ((void *) seg_info);
GMT_free ((void *) seg_area);
GMT_free ((void *) seg_start);
}
int GMT_init_br (char which, char res, struct GMT_BR *c, double w, double e, double s, double n)
/* which: r(iver) or b(order) */
/* res: Resolution (f, h, i, l, c */
{
int i, nb, idiv, iw, ie, is, in, this_south, this_west;
short *stmp;
int *itmp;
size_t start[1], count[1];
char file[32], path[BUFSIZ];
if (which == 'r')
sprintf (file, "binned_river_%c.cdf\0", res);
else
sprintf (file, "binned_border_%c.cdf\0", res);
if (!GMT_getpathname (file, path)) return (-1); /* Failed to find file */
check_nc_status (nc_open (path, NC_NOWRITE, &c->cdfid));
/* Get all id tags */
check_nc_status (nc_inq_varid (c->cdfid, "Bin_size_in_minutes", &c->bin_size_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_360_longitude_range", &c->bin_nx_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_180_degree_latitude_range", &c->bin_ny_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_bins_in_file", &c->n_bin_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_segments_in_file", &c->n_seg_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_points_in_file", &c->n_pt_id));
check_nc_status (nc_inq_varid (c->cdfid, "Id_of_first_segment_in_a_bin", &c->bin_firstseg_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_segments_in_a_bin", &c->bin_nseg_id));
check_nc_status (nc_inq_varid (c->cdfid, "N_points_for_a_segment", &c->seg_n_id));
check_nc_status (nc_inq_varid (c->cdfid, "Hierarchial_level_of_a_segment", &c->seg_level_id));
check_nc_status (nc_inq_varid (c->cdfid, "Id_of_first_point_in_a_segment", &c->seg_start_id));
check_nc_status (nc_inq_varid (c->cdfid, "Relative_longitude_from_SW_corner_of_bin", &c->pt_dx_id));
check_nc_status (nc_inq_varid (c->cdfid, "Relative_latitude_from_SW_corner_of_bin", &c->pt_dy_id));
/* Get attributes */
check_nc_status (nc_get_att_text (c->cdfid, c->pt_dx_id, "units", c->units));
check_nc_status (nc_get_att_text (c->cdfid, NC_GLOBAL, "title", c->title));
check_nc_status (nc_get_att_text (c->cdfid, NC_GLOBAL, "source", c->source));
/* Get global variables */
start[0] = 0;
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_size_id, start, &c->bin_size));
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_nx_id, start, &c->bin_nx));
check_nc_status (nc_get_var1_int (c->cdfid, c->bin_ny_id, start, &c->bin_ny));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_bin_id, start, &c->n_bin));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_seg_id, start, &c->n_seg));
check_nc_status (nc_get_var1_int (c->cdfid, c->n_pt_id, start, &c->n_pt));
c->scale = (c->bin_size / 60.0) / 65535.0;
c->bsize = c->bin_size / 60.0;
c->bins = (int *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (int), "GMT_init_br");
/* Round off area to nearest multiple of block-dimension */
iw = (int)(floor (w / c->bsize) * c->bsize);
ie = (int)(ceil (e / c->bsize) * c->bsize);
is = 90 - (int)(ceil ((90.0 - s) / c->bsize) * c->bsize);
in = 90 - (int)(floor ((90.0 - n) / c->bsize) * c->bsize);
idiv = irint (360.0 / c->bsize); /* Number of blocks per latitude band */
for (i = nb = 0; i < c->n_bin; i++) { /* Find which bins are needed */
this_south = 90 - (int)(c->bsize * ((i / idiv) + 1));
if (this_south < is || this_south >= in) continue;
this_west = (int)(c->bsize * (i % idiv)) - 360;
while (this_west < iw) this_west += 360;
if (this_west >= ie) continue;
c->bins[nb] = i;
nb++;
}
c->bins = (int *) GMT_memory ((void *)c->bins, (size_t)nb, sizeof (int), "GMT_init_br");
c->nb = nb;
/* Get bin variables, then extract only those corresponding to the bins to use */
/* Allocate space for arrays of bin information */
c->bin_nseg = (short *) GMT_memory (VNULL, (size_t)nb, sizeof (short), "GMT_init_br");
c->bin_firstseg = (int *) GMT_memory (VNULL, (size_t)nb, sizeof (int), "GMT_init_br");
count[0] = c->n_bin;
stmp = (short *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (short), "GMT_init_br");
check_nc_status (nc_get_vara_short (c->cdfid, c->bin_nseg_id, start, count, stmp));
for (i = 0; i < c->nb; i++) c->bin_nseg[i] = stmp[c->bins[i]];
GMT_free ((void *)stmp);
itmp = (int *) GMT_memory (VNULL, (size_t)c->n_bin, sizeof (int), "GMT_init_br");
check_nc_status (nc_get_vara_int (c->cdfid, c->bin_firstseg_id, start, count, itmp));
for (i = 0; i < c->nb; i++) c->bin_firstseg[i] = itmp[c->bins[i]];
GMT_free ((void *)itmp);
return (0);
}
void GMT_get_br_bin (int b, struct GMT_BR *c, int *level, int n_levels)
/* b: index number into c->bins */
/* level: Levels of features to extract */
/* n_levels: # of such levels. 0 means use all levels */
{
size_t start[1], count[1];
int *seg_start;
short *seg_n, *seg_level;
int s, i, k, skip;
c->lon_sw = (c->bins[b] % c->bin_nx) * c->bin_size / 60.0;
c->lat_sw = 90.0 - ((c->bins[b] / c->bin_nx) + 1) * c->bin_size / 60.0;
c->ns = c->bin_nseg[b];
if (c->ns == 0) return;
start[0] = c->bin_firstseg[b];
count[0] = c->bin_nseg[b];
seg_n = (short *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (short), "GMT_get_br_bin");
seg_level = (short *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (short), "GMT_get_br_bin");
seg_start = (int *) GMT_memory (VNULL, (size_t)c->bin_nseg[b], sizeof (int), "GMT_get_br_bin");
check_nc_status (nc_get_vara_short (c->cdfid, c->seg_n_id, start, count, seg_n));
check_nc_status (nc_get_vara_short (c->cdfid, c->seg_level_id, start, count, seg_level));
check_nc_status (nc_get_vara_int (c->cdfid, c->seg_start_id, start, count, seg_start));
c->seg = (struct GMT_BR_SEGMENT *) GMT_memory (VNULL, (size_t)c->ns, sizeof (struct GMT_BR_SEGMENT), "GMT_get_br_bin");
for (s = i = 0; i < c->ns; i++) {
if (n_levels == 0)
skip = FALSE;
else {
for (k = 0, skip = TRUE; skip && k < n_levels; k++)
if (seg_level[i] == level[k]) skip = FALSE;
}
if (skip) continue;
c->seg[s].n = seg_n[i];
c->seg[s].level = seg_level[i];
c->seg[s].dx = (short *) GMT_memory (VNULL, (size_t)c->seg[s].n, sizeof (short), "GMT_get_br_bin");
c->seg[s].dy = (short *) GMT_memory (VNULL, (size_t)c->seg[s].n, sizeof (short), "GMT_get_br_bin");
start[0] = seg_start[i];
count[0] = c->seg[s].n;
check_nc_status (nc_get_vara_short (c->cdfid, c->pt_dx_id, start, count, c->seg[s].dx));
check_nc_status (nc_get_vara_short (c->cdfid, c->pt_dy_id, start, count, c->seg[s].dy));
s++;
}
c->ns = s;
GMT_free ((void *) seg_n);
GMT_free ((void *) seg_level);
GMT_free ((void *) seg_start);
}
int GMT_assemble_shore (struct GMT_SHORE *c, int dir, int first_level, BOOLEAN assemble, BOOLEAN shift, double west, double east, struct POL **pol)
/* assemble: TRUE if polygons is needed */
/* shift: TRUE if longitudes may have to be shifted */
/* edge: Edge test for shifting */
{
struct POL *p;
int start_side, next_side, id, P = 0, more, p_alloc, wet_or_dry, use_this_level;
int n_alloc, cid, nid, add, first_pos, entry_pos, n, low_level, high_level;
BOOLEAN completely_inside;
double *xtmp, *ytmp, plon, plat;
if (!assemble) { /* Easy, just need to scale all segments to degrees and return */
p = (struct POL *) GMT_memory (VNULL, (size_t)c->ns, sizeof (struct POL), "GMT_assemble_shore");
for (id = 0; id < c->ns; id++) {
p[id].lon = (double *) GMT_memory (VNULL, (size_t)c->seg[id].n, sizeof (double), "GMT_assemble_shore");
p[id].lat = (double *) GMT_memory (VNULL, (size_t)c->seg[id].n, sizeof (double), "GMT_assemble_shore");
p[id].n = GMT_copy_to_shore_path (p[id].lon, p[id].lat, c, id);
p[id].level = c->seg[id].level;
p[id].interior = FALSE;
GMT_shore_path_shift2 (p[id].lon, p[id].lat, p[id].n, west, east, c->leftmost_bin);
}
*pol = p;
return (c->ns);
}
for (n = high_level = 0; n < 4; n++) high_level = MAX (c->node_level[n], high_level);
wet_or_dry = (dir == 1) ? 1 : 0;
use_this_level = (high_level%2 == wet_or_dry && high_level >= first_level);
if (c->ns == 0 && !use_this_level) return (0); /* No polygons for this bin */
/* Here we must assemble [at least one] polygons in the correct order */
for (n = 0, completely_inside = TRUE; completely_inside && n < c->ns; n++) if (c->seg[n].entry != 4) completely_inside = FALSE;
shore_prepare_sides (c, dir);
p_alloc = (c->ns == 0) ? 1 : GMT_SMALL_CHUNK;
p = (struct POL *) GMT_memory (VNULL, (size_t)p_alloc, sizeof (struct POL), "GMT_assemble_shore");
low_level = MAX_LEVEL;
if (completely_inside && use_this_level) { /* Must include path of this bin outline as first polygon */
p[0].n = GMT_graticule_path (&p[0].lon, &p[0].lat, dir, c->lon_corner[3], c->lon_corner[1], c->lat_corner[0], c->lat_corner[2]);
p[0].level = c->node_level[0]; /* Any corner will do */
p[0].interior = FALSE;
P = 1;
}
while (c->n_entries > 0) { /* More segments to connect */
low_level = MAX_LEVEL;
start_side = 0;
id = GMT_shore_get_first_entry (c, dir, &start_side);
next_side = c->seg[id].exit;
n_alloc = c->seg[id].n;
p[P].lon = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].lat = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
n = GMT_copy_to_shore_path (p[P].lon, p[P].lat, c, id);
if ((int)c->seg[id].level < low_level) low_level = (int)c->seg[id].level;
more = TRUE;
first_pos = GMT_shore_get_position (start_side, c->seg[id].dx[0], c->seg[id].dy[0]);
while (more) {
id = GMT_shore_get_next_entry (c, dir, next_side, id);
if (id < 0) { /* Corner */
cid = id + 4;
nid = (dir == 1) ? (cid + 1) % 4 : cid;
if ((add = GMT_map_path (p[P].lon[n-1], p[P].lat[n-1], c->lon_corner[cid], c->lat_corner[cid], &xtmp, &ytmp))) {
n_alloc += add;
p[P].lon = (double *) GMT_memory ((void *)p[P].lon, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].lat = (double *) GMT_memory ((void *)p[P].lat, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
memcpy ((void *)&p[P].lon[n], (void *)xtmp, (size_t)(add * sizeof (double)));
memcpy ((void *)&p[P].lat[n], (void *)ytmp, (size_t)(add * sizeof (double)));
n += add;
}
next_side = ((id + 4) + dir + 4) % 4;
if ((int)c->node_level[nid] < low_level) low_level = (int)c->node_level[nid];
}
else {
GMT_shore_to_degree (c, c->seg[id].dx[0], c->seg[id].dy[0], &plon, &plat);
if ((add = GMT_map_path (p[P].lon[n-1], p[P].lat[n-1], plon, plat, &xtmp, &ytmp))) {
n_alloc += add;
p[P].lon = (double *) GMT_memory ((void *)p[P].lon, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].lat = (double *) GMT_memory ((void *)p[P].lat, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
memcpy ((void *)&p[P].lon[n], (void *)xtmp, (size_t)(add * sizeof (double)));
memcpy ((void *)&p[P].lat[n], (void *)ytmp, (size_t)(add * sizeof (double)));
n += add;
}
entry_pos = GMT_shore_get_position (next_side, c->seg[id].dx[0], c->seg[id].dy[0]);
if (next_side == start_side && entry_pos == first_pos)
more = FALSE;
else {
n_alloc += c->seg[id].n;
p[P].lon = (double *) GMT_memory ((void *)p[P].lon, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].lat = (double *) GMT_memory ((void *)p[P].lat, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
n += GMT_copy_to_shore_path (&p[P].lon[n], &p[P].lat[n], c, id);
next_side = c->seg[id].exit;
if ((int)c->seg[id].level < low_level) low_level = (int)c->seg[id].level;
}
}
if (add) {
GMT_free ((void *)xtmp);
GMT_free ((void *)ytmp);
}
}
p[P].n = n;
p[P].interior = FALSE;
p[P].level = (dir == 1) ? 2 * ((low_level - 1) / 2) + 1: 2 * (low_level/2);
P++;
if (P == p_alloc) {
p_alloc += GMT_SMALL_CHUNK;
p = (struct POL *) GMT_memory ((void *)p, (size_t)p_alloc, sizeof (struct POL), "GMT_assemble_shore");
}
}
/* Then add all interior polygons, if any */
for (id = 0; id < c->ns; id++) {
if (c->seg[id].entry < 4) continue;
n_alloc = c->seg[id].n;
p[P].lon = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].lat = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), "GMT_assemble_shore");
p[P].n = GMT_copy_to_shore_path (p[P].lon, p[P].lat, c, id);
p[P].interior = TRUE;
p[P].level = c->seg[id].level;
P++;
if (P == p_alloc) {
p_alloc += GMT_SMALL_CHUNK;
p = (struct POL *) GMT_memory ((void *)p, (size_t)p_alloc, sizeof (struct POL), "GMT_assemble_shore");
}
}
GMT_shore_pau_sides (c);
if (c->ns > 0) p = (struct POL *) GMT_memory ((void *)p, (size_t)P, sizeof (struct POL), "GMT_assemble_shore");
for (id = 0; id < P; id++) GMT_shore_path_shift2 (p[id].lon, p[id].lat, p[id].n, west, east, c->leftmost_bin);
*pol = p;
return (P);
}
int GMT_assemble_br (struct GMT_BR *c, BOOLEAN shift, double edge, struct POL **pol)
/* shift: TRUE if longitudes may have to be shifted */
/* edge: Edge test for shifting */
{
struct POL *p;
int id;
p = (struct POL *) GMT_memory (VNULL, (size_t)c->ns, sizeof (struct POL), "GMT_assemble_br");
for (id = 0; id < c->ns; id++) {
p[id].lon = (double *) GMT_memory (VNULL, (size_t)c->seg[id].n, sizeof (double), "GMT_assemble_br");
p[id].lat = (double *) GMT_memory (VNULL, (size_t)c->seg[id].n, sizeof (double), "GMT_assemble_br");
p[id].n = GMT_copy_to_br_path (p[id].lon, p[id].lat, c, id);
p[id].level = c->seg[id].level;
if (shift) GMT_shore_path_shift (p[id].lon, p[id].lat, p[id].n, edge);
}
*pol = p;
return (c->ns);
}
void GMT_free_shore (struct GMT_SHORE *c)
{ /* Removes allocated variables for this block only */
int i;
for (i = 0; i < c->ns; i++) {
GMT_free ((void *)c->seg[i].dx);
GMT_free ((void *)c->seg[i].dy);
}
if (c->ns) GMT_free ((void *)c->seg);
}
void GMT_free_br (struct GMT_BR *c)
{ /* Removes allocated variables for this block only */
int i;
for (i = 0; i < c->ns; i++) {
GMT_free ((void *)c->seg[i].dx);
GMT_free ((void *)c->seg[i].dy);
}
if (c->ns) GMT_free ((void *)c->seg);
}
void GMT_shore_cleanup (struct GMT_SHORE *c)
{
GMT_free ((void *)c->bins);
GMT_free ((void *)c->bin_info);
GMT_free ((void *)c->bin_nseg);
GMT_free ((void *)c->bin_firstseg);
check_nc_status (nc_close (c->cdfid));
}
void GMT_br_cleanup (struct GMT_BR *c)
{
GMT_free ((void *)c->bins);
GMT_free ((void *)c->bin_nseg);
GMT_free ((void *)c->bin_firstseg);
check_nc_status (nc_close (c->cdfid));
}
/* ---------- LOWER LEVEL FUNCTIONS CALLED BY THE ABOVE ------------ */
int GMT_copy_to_shore_path (double *lon, double *lat, struct GMT_SHORE *s, int id)
{
int i;
for (i = 0; i < (int)s->seg[id].n; i++)
GMT_shore_to_degree (s, s->seg[id].dx[i], s->seg[id].dy[i], &lon[i], &lat[i]);
return (s->seg[id].n);
}
int GMT_copy_to_br_path (double *lon, double *lat, struct GMT_BR *s, int id)
{
int i;
for (i = 0; i < (int)s->seg[id].n; i++)
GMT_br_to_degree (s, s->seg[id].dx[i], s->seg[id].dy[i], &lon[i], &lat[i]);
return (s->seg[id].n);
}
void GMT_shore_to_degree (struct GMT_SHORE *c, short int dx, short int dy, double *lon, double *lat)
{
*lon = c->lon_sw + ((unsigned short)dx) * c->scale;
*lat = c->lat_sw + ((unsigned short)dy) * c->scale;
}
void GMT_br_to_degree (struct GMT_BR *c, short int dx, short int dy, double *lon, double *lat)
{
*lon = c->lon_sw + ((unsigned short)dx) * c->scale;
*lat = c->lat_sw + ((unsigned short)dy) * c->scale;
}
int GMT_shore_get_next_entry (struct GMT_SHORE *c, int dir, int side, int id)
{
/* Finds the next entry point on the given side that is further away
* in the <dir> direction than previous point. It removes the info
* regarding the new entry from the SIDE structure */
int k, pos, n;
if (id < 0)
pos = (dir == 1) ? 0 : MAX_DELTA;
else {
n = c->seg[id].n - 1;
pos = GMT_shore_get_position (side, c->seg[id].dx[n], c->seg[id].dy[n]);
}
if (dir == 1) {
for (k = 0; k < (int)c->nside[side] && (int)c->side[side][k].pos < pos; k++);
id = c->side[side][k].id;
for (k++; k < c->nside[side]; k++) c->side[side][k-1] = c->side[side][k];
c->nside[side]--;
}
else {
for (k = 0; k < (int)c->nside[side] && (int)c->side[side][k].pos > pos; k++);
id = c->side[side][k].id;
for (k++; k < c->nside[side]; k++) c->side[side][k-1] = c->side[side][k];
c->nside[side]--;
}
if (id >= 0) c->n_entries--;
return (id);
}
int GMT_shore_get_first_entry (struct GMT_SHORE *c, int dir, int *side)
{
int try = 0;
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 */
try++;
*side = (*side + dir + 4) % 4;
}
if (try == 4) return (-5);
return (c->side[*side][0].id);
}
int GMT_shore_asc_sort (const void *a, const void *b)
{
if (((struct SIDE *)a)->pos < ((struct SIDE *)b)->pos) return (-1);
if (((struct SIDE *)a)->pos > ((struct SIDE *)b)->pos) return (1);
return (0);
}
int GMT_shore_desc_sort (const void *a, const void *b)
{
if (((struct SIDE *)a)->pos < ((struct SIDE *)b)->pos) return (1);
if (((struct SIDE *)a)->pos > ((struct SIDE *)b)->pos) return (-1);
return (0);
}
void GMT_shore_pau_sides (struct GMT_SHORE *c)
{
int i;
for (i = 0; i < 4; i++) GMT_free ((void *)c->side[i]);
}
void GMT_free_polygons (struct POL *p, int n)
{
int k;
for (k = 0; k < n; k++) {
GMT_free ((void *)p[k].lon);
GMT_free ((void *)p[k].lat);
}
}
void GMT_shore_path_shift (double *lon, double *lat, int n, double edge)
{
int i;
for (i = 0; i < n; i++) if (lon[i] >= edge) lon[i] -= 360.0;
}
void GMT_shore_path_shift2old (double *lon, double *lat, int n, double west, double east)
{
int i;
/* for (i = 0; i < n; i++) if (lon[i] >= east && (lon[i]-360.0) > west) lon[i] -= 360.0; */
for (i = 0; i < n; i++) if (lon[i] > east && (lon[i]-360.0) >= west) lon[i] -= 360.0;
}
void GMT_shore_path_shift2 (double *lon, double *lat, int n, double west, double east, BOOLEAN leftmost)
{
int i;
if (leftmost) { /* Must check this bin differently */
for (i = 0; i < n; i++) if (lon[i] >= east && (lon[i]-360.0) >= west) lon[i] -= 360.0;
}
else {
for (i = 0; i < n; i++) if (lon[i] > east && (lon[i]-360.0) >= west) lon[i] -= 360.0;
}
}
int GMT_shore_get_position (int side, short int x, short int y)
{
/* Returns the position along the given side */
return ((side%2) ? ((side == 1) ? (unsigned short)y : MAX_DELTA - (unsigned short)y) : ((side == 0) ? (unsigned short)x : MAX_DELTA - (unsigned short)x));
}
void shore_prepare_sides (struct GMT_SHORE *c, int dir)
{
int s, i, n[4];
c->lon_corner[0] = c->lon_sw + ((dir == 1) ? c->bsize : 0.0);
c->lon_corner[1] = c->lon_sw + c->bsize;
c->lon_corner[2] = c->lon_sw + ((dir == 1) ? 0.0 : c->bsize);
c->lon_corner[3] = c->lon_sw;
c->lat_corner[0] = c->lat_sw;
c->lat_corner[1] = c->lat_sw + ((dir == 1) ? c->bsize : 0.0);
c->lat_corner[2] = c->lat_sw + c->bsize;
c->lat_corner[3] = c->lat_sw + ((dir == 1) ? 0.0 : c->bsize);
for (i = 0; i < 4; i++) c->nside[i] = n[i] = 1;
/* for (s = 0; s < c->ns; s++) if (c->seg[s].level < 3 && c->seg[s].entry < 4) c->nside[c->seg[s].entry]++; */
for (s = 0; s < c->ns; s++) if (c->seg[s].entry < 4) c->nside[c->seg[s].entry]++;
for (i = c->n_entries = 0; i < 4; i++) { /* Allocate memory and add corners */
c->side[i] = (struct SIDE *) GMT_memory (VNULL, (size_t)c->nside[i], sizeof (struct SIDE), "shore_prepare_sides");
c->side[i][0].pos = (dir == 1) ? MAX_DELTA : 0;
c->side[i][0].id = i - 4;
c->n_entries += c->nside[i] - 1;
}
for (s = 0; s < c->ns; s++) { /* Add entry points */
/* if (c->seg[s].level > 2 || (i = c->seg[s].entry) == 4) continue; */
if ((i = c->seg[s].entry) == 4) continue;
c->side[i][n[i]].pos = GMT_shore_get_position (i, c->seg[s].dx[0], c->seg[s].dy[0]);
c->side[i][n[i]].id = s;
n[i]++;
}
for (i = 0; i < 4; i++) { /* sort on position */
if (dir == 1)
qsort ((void *)c->side[i], (size_t)c->nside[i], sizeof (struct SIDE), GMT_shore_asc_sort);
else
qsort ((void *)c->side[i], (size_t)c->nside[i], sizeof (struct SIDE), GMT_shore_desc_sort);
}
}