Address issue #159 (#163)

We never bothered to clip the output of -M to fit the given -R so data would exceed those limits. THis was partly because it did not seem to be a big deal, but it isnt very nice. Now we clip.
GenericMappingTools · Nov 5, 2018 · c5a5435 · c5a5435
1 parent 69d48ad
commit c5a5435
Show file tree

Hide file tree

Showing 3 changed files with 1,292 additions and 1,853 deletions.
diff --git a/src/gmt_map.c b/src/gmt_map.c
@@ -1663,7 +1663,7 @@ GMT_LOCAL uint64_t map_wesn_clip_old (struct GMT_CTRL *GMT, double *lon, double
 }
 
 /*! . */
-GMT_LOCAL uint64_t map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t n_orig, double **x, double **y, uint64_t *total_nx) {
+uint64_t map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t n_orig, double **x, double **y, uint64_t *total_nx) {
 	char *x_type = NULL;
 	size_t n_alloc = 0, n_x_alloc = 0, n_t_alloc = 0;
 	uint64_t new_n, i, n_get, n, m, n_cross = 0, *x_index = NULL;

diff --git a/src/pscoast.c b/src/pscoast.c
@@ -691,6 +691,8 @@ int GMT_coast (void *V_API, int mode, void *args) {
 	return GMT_pscoast (V_API, mode, args);
 }
 
+EXTERN_MSC uint64_t map_wesn_clip (struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t n_orig, double **x, double **y, uint64_t *total_nx);
+
 int GMT_pscoast (void *V_API, int mode, void *args) {
 	/* High-level function that implements the pscoast task */
 
@@ -1074,7 +1076,9 @@ int GMT_pscoast (void *V_API, int mode, void *args) {
 			}
 
 			for (i = 0; i < np; i++) {
-				if (Ctrl->M.active) {
+				if (Ctrl->M.active) {	/* Clip to specified region - this gives x,y coordinates in inches */
+					uint64_t unused, n_out = map_wesn_clip (GMT, p[i].lon, p[i].lat, p[i].n, &xtmp, &ytmp, &unused);
+
 					if (!Ctrl->M.single) {
 						sprintf (GMT->current.io.segment_header, "Shore Bin # %d, Level %d", bin, p[i].level);
 						GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
@@ -1083,11 +1087,12 @@ int GMT_pscoast (void *V_API, int mode, void *args) {
 						out[GMT_X] = out[GMT_Y] = GMT->session.d_NaN;
 						GMT_Put_Record (API, GMT_WRITE_DATA, Out);
 					}
-					for (k = 0; k < p[i].n; k++) {
-						out[GMT_X] = p[i].lon[k];
-						out[GMT_Y] = p[i].lat[k];
+					for (k = 0; k < (int)n_out; k++) {	/* Convert back to lon,lat before writing output */
+						gmt_xy_to_geo (GMT, &out[GMT_X], &out[GMT_Y], xtmp[k], ytmp[k]);
 						GMT_Put_Record (API, GMT_WRITE_DATA, Out);
 					}
+					gmt_M_free (GMT, xtmp);
+					gmt_M_free (GMT, ytmp);
 				}
 				else if (Ctrl->W.use[p[i].level-1]) {
 					if (donut_hell) p[i].n = (int)gmt_fix_up_path (GMT, &p[i].lon, &p[i].lat, p[i].n, 0.0, 0);