Function gmtplot_geo_polygon_segment

src/gmt_plot.c:4722–4797

}

GMT_LOCAL uint64_t gmtplot_geo_polygon_segment (struct GMT_CTRL *GMT, struct GMT_DATASEGMENT *S, bool add_pole, bool first, const char *comment) {
	/* Handles the laying down of polygons suitable for filling only; outlines are done separately later.
	 * Polar caps need special treatment in that we must add a detour to the pole.
	 * That detour will not be drawn, only used for fill. However, due to the insanity that is called GIS,
	 * some user polygons may already have artificial lines drawn to the pole in order to work in GIS.
	 * Thus, if we detect such a pole as part of the line then we do NOT add yet another detour.  */

	uint64_t n = S->n_rows, k;
	double *plon = S->data[GMT_X], *plat = S->data[GMT_Y], t_lat;	/* Default is to plot incoming array as is via plon,plat pointers */
	bool ap = gmtplot_at_pole (plat, n);	/* Is the first and last point exactly at the pole? */
	bool free_memory = false;
	struct GMT_DATASEGMENT_HIDDEN *SH = gmt_get_DS_hidden (S);
	if (ap) plon[n-1] = plon[0];	/* Just enforce the same longitude at the pole point */
	GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Polar cap: %d\n", (int)add_pole);
	if (add_pole) {	/* Make sure there is not already a detour in the data as given */
		double p_lat = SH->pole * 90.0;	/* Latitude of the pole in question */
		bool need_detour = true;	/* Until proven otherwise we assume we must add a detour */
		if (GMT->common.R.oblique) {	/* Determine if any of the 4 map corners are inside this polygon */
			double X, Y;
			gmt_set_inside_mode (GMT, NULL, GMT_IOO_SPHERICAL);
			gmt_xy_to_geo (GMT, &X, &Y, GMT->current.proj.rect[XLO], GMT->current.proj.rect[YLO]);
			GMT->current.proj.corner[0] = gmt_inonout (GMT, X, Y, S);
			gmt_xy_to_geo (GMT, &X, &Y, GMT->current.proj.rect[XHI], GMT->current.proj.rect[YLO]);
			GMT->current.proj.corner[1] = gmt_inonout (GMT, X, Y, S);
			gmt_xy_to_geo (GMT, &X, &Y, GMT->current.proj.rect[XHI], GMT->current.proj.rect[YHI]);
			GMT->current.proj.corner[2] = gmt_inonout (GMT, X, Y, S);
			gmt_xy_to_geo (GMT, &X, &Y, GMT->current.proj.rect[XLO], GMT->current.proj.rect[YHI]);
			GMT->current.proj.corner[3] = gmt_inonout (GMT, X, Y, S);
			need_detour = false;	/* Trying a different tack for these cases */
		}
		for (k = 0; need_detour && k < S->n_rows; k++) {	/* Check every point */
			if (doubleAlmostEqual (S->data[GMT_Y][k], p_lat)) {	/* Point is exactly at the pole in question */
				/* We want to distinguish between a path that gently touches the pole and one that has a fake straight detour to the pole.
				 * We assume a fake detour will have the same longitudes for this point and the previous and that they are both either +/-180 or 0. */
				if (k && doubleAlmostEqual (S->data[GMT_X][k], S->data[GMT_X][k-1]) && (doubleAlmostEqual (fabs (S->data[GMT_X][k]), 180.0) || gmt_M_is_zero (S->data[GMT_X][k])))
					need_detour = false;	/* Well, what do you know. Probably arcGIS or some other handicapped program */
			}
		}
		if (!need_detour) {	/* Do not add another detour but process via gmt_geo_polarcap_segment to handle jumps in our polygon */
			GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Path already had a detour to the pole, skip adding another detour\n");
			add_pole = false;
			n = gmt_geo_polarcap_segment (GMT, S, &plon, &plat);
			if (plon != S->data[GMT_X])	/* Pointer changed so we allocated, must delete below */
				free_memory = true;
		}
	}
	if (add_pole) {	/* If we get here then a detour will be needed */
		if ((n = gmt_geo_polarcap_segment (GMT, S, &plon, &plat)) == 0) {	/* Not a global map */
			/* Here we must detour to the N or S pole, then resample the path */
			n = S->n_rows + 2;	/* Add new first and last point to connect to the pole */
			plon = gmt_M_memory (GMT, NULL, n, double);	/* This memory must be freed below */
			plat = gmt_M_memory (GMT, NULL, n, double);
			free_memory = true;
			t_lat = SH->pole * 90.0;	/* This is presumably the correct pole, but could fail if just touching the pole */
			if (S->data[GMT_Y][0] * t_lat < 0.0) t_lat = -t_lat;	/* Well, I'll be damned... */
			plat[0] = plat[n-1] = t_lat;
			plon[0] = S->data[GMT_X][0];
			plon[n-1] = S->data[GMT_X][S->n_rows-1];