Final version of updates for 6.5: grdimage, colorbar, grdmix (#8251)

* Update grdimage.c

Moving parts over from old branch

* More brought over and still no failures

* Update grdimage.c

* Update grdimage.c

* More copying with no failures

* More OK updates

* Update grdimage.c

* Add test grid_transp.sh

* Handle -Q+i for image

* Add transp_mix test and PS

* Update grdimage.c

Did not test OpenMP so missed some new variables.

* Update transp_mix.sh

* Update grdimage.dvc

* Fix OpenMP issue

* Add image_vartrans.sh test and PS to DVC

* Update grdimage.c

* Update image_vartrans.sh

* polish

* Add test grdimage_Q_effects.sh and PS to DVC

* indents

* Update grdimage.c

* Be more careful on memory layout for iimages

This PR does now write a blank space after TRB, but places a there only if the images has transparency.

* Fill in missing variables in OPENMP pragmas

Especially in grdmix

* Add gmt_setrgb and remove OPENMP loop when PSL memory is required

* Update grdimage.c

* Update grdimage.c

* Fix opacity issue

* Update psscale.c

Finalize the -Li|I positioning and frame dimensions and not break previous tests.

* Add new test and in DVC

* alpha_count instead of alpha, which is unclear.

---------

Co-authored-by: Joaquim <jluis@ualg.pt>

src/gmt_hidden.h

@@ -177,6 +177,7 @@ struct GMT_GRID_HEADER_HIDDEN {
 	char flags[4];                        /* Flags used for ESRI grids */
 	char *pocket;                         /* GDAL: A working variable handy to transmit info between funcs e.g. +b<band_info> to gdalread */
 	double bcr_threshold;                 /* sum of cardinals must >= threshold in bilinear; else NaN */
+	unsigned int used_indexed_rgb;        /* 1 if originally an indexed RGB, 0 otherwise */
 	unsigned int has_NaN_rgb;             /* Is 1 if an indexed RGB image had a specific NaN color for transparency */
 	unsigned int has_NaNs;                /* Is 2 if the grid contains any NaNs, 1 if it does not, and 0 if no check has yet happened */
 	unsigned int bcr_interpolant;         /* Interpolation function used (0, 1, 2, 3) */

src/gmt_plot.c

@@ -6914,6 +6914,12 @@ void gmt_BB_clip_on (struct GMT_CTRL *GMT, double rgb[], unsigned int flag) {
 	PSL_beginclipping (PSL, work_x, work_y, 5, rgb, flag);
 }
 
+void gmt_setrgb (struct GMT_CTRL *GMT, double *rgb) {
+	struct PSL_CTRL *PSL= GMT->PSL;
+	/* Fill with a color */
+	PSL_setrgb (PSL, rgb);
+}
+
 void gmt_setfill (struct GMT_CTRL *GMT, struct GMT_FILL *fill, int outline) {
 	struct PSL_CTRL *PSL= GMT->PSL;
 	if (!fill) /* NO fill pointer = no fill */

src/gmt_prototypes.h

@@ -304,6 +304,7 @@ EXTERN_MSC void gmt_BB_clip_on (struct GMT_CTRL *GMT, double rgb[], unsigned int
 EXTERN_MSC void gmt_plot_line (struct GMT_CTRL *GMT, double *x, double *y, unsigned int *pen, uint64_t n, unsigned int mode);
 EXTERN_MSC void gmt_setpen (struct GMT_CTRL *GMT, struct GMT_PEN *pen);
 EXTERN_MSC void gmt_setfill (struct GMT_CTRL *GMT, struct GMT_FILL *fill, int outline);
+EXTERN_MSC void gmt_setrgb (struct GMT_CTRL *GMT, double *rgb);
 EXTERN_MSC void gmt_vertical_axis (struct GMT_CTRL *GMT, unsigned int mode);
 EXTERN_MSC void gmt_xy_axis (struct GMT_CTRL *GMT, double x0, double y0, double length, double val0, double val1, struct GMT_PLOT_AXIS *A, bool below, unsigned int side);
 EXTERN_MSC void gmt_xy_axis2 (struct GMT_CTRL *GMT, double x0, double y0, double length, double val0, double val1, struct GMT_PLOT_AXIS *A, bool below, bool annotate, unsigned side);

src/grdmix.c

@@ -47,7 +47,8 @@
 
 struct GRDMIX_AIW {	/* For various grid, image, or constant arguments */
 	bool active;
-	bool opacity;	/* true if we hvae opacity instead of transparency [Default] */
+	/* Our implementation is a bit backwards, so to get expected results we turn opacity on by default and let +o reverse that */
+	bool opacity;	/* true if we have opacity instead of transparency [Default] */
 	unsigned int mode;	/* 0 a file, 1 a constant */
 	char *file;
 	double value;
@@ -580,7 +581,7 @@ EXTERN_MSC int GMT_grdmix (void *V_API, int mode, void *args) {
 			}
 			else {
 #ifdef _OPENMP
-#pragma omp parallel for private(row,col,node) shared(GMT,G,off,scale,I_in)
+#pragma omp parallel for private(row,col,node) shared(GMT,G,off,I_in)
 #endif
 				gmt_M_grd_loop (GMT, G, row, col, node)
 					G->data[node] = I_in[0]->data[node+off];
@@ -627,9 +628,9 @@ EXTERN_MSC int GMT_grdmix (void *V_API, int mode, void *args) {
 					GMT_Report (API, GMT_MSG_WARNING, "Component grid values in %s exceed 0-1 range, probably need to specify -Ni\n", Ctrl->In.file[band]);
 			}
 			if (Ctrl->I.active && Ctrl->In.n_in == 3) {	/* Make the most work-intensive version under OpenMP */
-	#ifdef _OPENMP
-	#pragma omp parallel for private(row,col,node,band,rgb,pix) shared(GMT,I,G_in,H,intens)
-	#endif
+#ifdef _OPENMP
+#pragma omp parallel for private(row,col,node,band,rgb,pix) shared(GMT,I,P,G_in,H,intens)
+#endif
 				gmt_M_grd_loop (GMT, I, row, col, node) {	/* The node is one per pixel in a band, so stride into additional bands */
 					if (P)	/* Get r/g/b from grid z-value via CPT lookup */
 						(void)gmt_get_rgb_from_z (GMT, P, G_in[0]->data[node], rgb);
@@ -769,11 +770,11 @@ EXTERN_MSC int GMT_grdmix (void *V_API, int mode, void *args) {
 			Return (GMT_RUNTIME_ERROR);
 		}
 #ifdef _OPENMP
-#pragma omp parallel for private(node) shared(H,I,alpha)
+#pragma omp parallel for private(node,transparency) shared(H,I,alpha,Ctrl)
 #endif
 		for (node = 0; node < (int64_t)H->size; node++)	{	/* Scale to 0-255 range */
 			transparency = gmt_M_is_dnan (alpha[node]) ? 1.0 : alpha[node];	/* NaN means full transparency */
-			if (Ctrl->A.opacity) transparency = 1.0 - alpha[node];	/* Turns out we got opacities */
+			if (!Ctrl->A.opacity) transparency = 1.0 - alpha[node];	/* Turns out we got opacities */
 			I->alpha[node] = gmt_M_u255 (transparency);
 		}
 		gmt_M_free (GMT, alpha);

src/postscriptlight.c

@@ -4982,6 +4982,22 @@ int PSL_setcolor (struct PSL_CTRL *PSL, double rgb[], int mode) {
 	return (PSL_NO_ERROR);
 }
 
+int PSL_setrgb (struct PSL_CTRL *PSL, double rgb[]) {
+	/* Set the fill (PSL_IS_FILL) r/g/b color
+	 * rgb[0] >= 0: rgb is the color with R G B in 0-1 range.
+	 */
+	if (!rgb) return (PSL_NO_ERROR);	/* NULL args to be ignored */
+	if (PSL_same_rgb (rgb, PSL->current.rgb[PSL_IS_FILL])) return (PSL_NO_ERROR);	/* Same color as already set */
+
+	/* Then, finally, set the color using psl_putcolor */
+	PSL_command (PSL, "{%s} FS\n", psl_putcolor (PSL, rgb, 0));
+
+	/* Update the current stroke/fill color information */
+
+	PSL_rgb_copy (PSL->current.rgb[PSL_IS_FILL], rgb);
+	return (PSL_NO_ERROR);
+}
+
 char * PSL_makepen (struct PSL_CTRL *PSL, double linewidth, double rgb[], char *pattern, double offset) {
 	/* Creates a text string with the corresponding PS command to set the pen */
 	static char buffer[PSL_BUFSIZ];

src/postscriptlight.h

@@ -490,6 +490,7 @@ EXTERN_MSC int PSL_setdefaults (struct PSL_CTRL *PSL, double xyscales[], double
 EXTERN_MSC int PSL_setdash (struct PSL_CTRL *PSL, char *pattern, double offset);
 EXTERN_MSC int PSL_setfill (struct PSL_CTRL *PSL, double rgb[], int outline);
 EXTERN_MSC int PSL_setfont (struct PSL_CTRL *PSL, int font_no);
+EXTERN_MSC int PSL_setrgb (struct PSL_CTRL *PSL, double rgb[]);
 EXTERN_MSC int PSL_setfontdims (struct PSL_CTRL *PSL, double supsub, double scaps, double sup_lc, double sup_uc, double sdown);
 EXTERN_MSC int PSL_setformat (struct PSL_CTRL *PSL, int n_decimals);
 EXTERN_MSC int PSL_setimage (struct PSL_CTRL *PSL, int image_no, char *imagefile, unsigned char *image, int image_dpi, unsigned int dim[], double f_rgb[], double b_rgb[]);

src/psscale.c

@@ -184,7 +184,6 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 1, "\n-D%s[+w<length>[/<width>]][+e[b|f][<length>]][+h|v][+j<justify>][+m[a|c|l|u]][+n[<txt>]]%s[+r]]", GMT_XYANCHOR, GMT_OFFSET);
 	GMT_Usage (API, -2, "Specify position and dimensions of the scale bar [JBC]. ");
 	gmt_refpoint_syntax (API->GMT, "D", NULL, GMT_ANCHOR_COLORBAR, 3);
-	//gmt_refpoint_syntax (API->GMT, "  ", "  Specify position and dimensions of the scale bar [JBC].", GMT_ANCHOR_COLORBAR, 1);
 	GMT_Usage (API, -2, "For -DJ|j w/TC|BC|ML|MR the values for +w (%d%% of map width), +h|v,+j,+o,+m have defaults. "
 		"You can override any of these settings with these explicit modifiers:", PSSCALE_L_SCALE);
 	GMT_Usage (API, 3, "+h Select a horizontal scale.");
@@ -1092,7 +1091,7 @@ GMT_LOCAL void psscale_draw_colorbar (struct GMT_CTRL *GMT, struct PSSCALE_CTRL
 			annot_off += bar_tick_len;
 			/* Extend x clearance by annotation width */
 			annot_off += hor_annot_width;
-			if (Ctrl->L.interval) annot_off += 0.4 * hor_annot_width;
+			if (Ctrl->L.interval) annot_off += 0.50 * hor_annot_width;
 			/* Increase width if there is a label */
 			if (label[0])
 				label_off = MAX (0.0, GMT->current.setting.map_label_offset[GMT_Y]) + GMT->current.setting.font_label.size / PSL_POINTS_PER_INCH;
@@ -1584,6 +1583,13 @@ GMT_LOCAL void psscale_draw_colorbar (struct GMT_CTRL *GMT, struct PSSCALE_CTRL
 		gmt_M_memset (text, 256U, char);
 		gmt_M_memset (test, 256U, char);
 		if (center && Ctrl->L.interval) {
+			sprintf (text, "%ld - %ld", lrint (floor (P->data[0].z_low)), lrint (ceil (P->data[0].z_high)));
+			sprintf (test, "%ld - %ld", lrint (floor (P->data[P->n_colors-1].z_low)), lrint (ceil (P->data[P->n_colors-1].z_high)));
+			hor_annot_width = ((MAX ((int)strlen (text), (int)strlen (test)) + 2*ndec) * GMT_DEC_WIDTH - 0.4 +
+				((ndec > 0) ? 2*GMT_PER_WIDTH : 0.0))
+				* GMT->current.setting.font_annot[GMT_PRIMARY].size * GMT->session.u2u[GMT_PT][GMT_INCH];
+		}
+		else if (center && Ctrl->L.interval && Ctrl->F.active) {
 			sprintf (format2, "%s%c%s", one_format, endash, one_format);
 			sprintf (text, format2, P->data[0].z_low, P->data[0].z_high);
 			sprintf (test, format2, P->data[P->n_colors-1].z_low, P->data[P->n_colors-1].z_high);
@@ -1606,8 +1612,8 @@ GMT_LOCAL void psscale_draw_colorbar (struct GMT_CTRL *GMT, struct PSSCALE_CTRL
 			y_annot = y_base + dir * (((len > 0.0) ? len : 0.0) + GMT->current.setting.map_annot_offset[GMT_PRIMARY] * cosd (Ctrl->S.angle));
 			justify = l_justify = (dir == -1) ? PSL_ML : PSL_MR;
 		}
-		else if (Ctrl->L.interval)
-			y_annot = y_base + dir * annot_off * 0.65;
+		else if (Ctrl->L.interval && Ctrl->F.active)
+			y_annot = y_base + dir * annot_off * 0.9;
 		else
 			y_annot = y_base + dir * annot_off;
 		if ((flip & PSSCALE_FLIP_ANNOT) == (flip & PSSCALE_FLIP_LABEL) / 2) y_label = y_base + dir * label_off;
@@ -1773,6 +1779,7 @@ GMT_LOCAL void psscale_draw_colorbar (struct GMT_CTRL *GMT, struct PSSCALE_CTRL
 					if (use_labels && (no_B_mode & PSSCALE_ANNOT_CUSTOM))
 						this_just = psscale_set_custom_annot (GMT, P, i, justify, l_justify, text);
 					else if (center && Ctrl->L.interval) {
+						sprintf (text, format, P->data[i].z_low, P->data[i].z_high);
 						if (Ctrl->L.interval == 2 && i == 0) {
 							sprintf (format2, "< %s", one_format);
 							sprintf (text, format2, P->data[i].z_high);

test/baseline/grdimage.dvc

@@ -1,5 +1,5 @@
 outs:
-- md5: d1c15fadf2e09ad118cc77c53bbbff92.dir
-  nfiles: 31
+- md5: 226f00f8bb56bfda27786069c3ecfa19.dir
+  nfiles: 35
   path: grdimage
   hash: md5

test/baseline/grdmix.dvc

@@ -1,6 +1,5 @@
 outs:
-- md5: 4aa0bfc2018d514b2b1e7af6231f457f.dir
-  size: 2345366
+- md5: 6549edda1e5c0c5ac5f3917baf574a83.dir
   nfiles: 2
   path: grdmix
   hash: md5

test/baseline/psscale.dvc

@@ -1,5 +1,6 @@
 outs:
-- md5: 452ffd60491b33499fed19555a0bb535.dir
-  size: 1400714
-  nfiles: 25
+- md5: 6e77da8a8cb4a2816e41162fe275360d.dir
+  size: 1431120
+  nfiles: 26
   path: psscale
+  hash: md5

test/grdimage/grdimage_Q_effects.sh

@@ -0,0 +1,35 @@
+#!/usr/bin/env bash
+#
+# Testing a simple grid with -Q options for transparency based on
+# a computed color (via grid z-value and CPT) or a specific color.
+
+gmt begin grdimage_Q_effects
+	gmt set FONT_TAG 12p PS_MEDIA letter
+	gmt grdmath -R0/100/0/50 -I1 -rp X 10 DIV FLOOR = stripes.nc
+	gmt grdmath -Rstripes.nc X 100 DIV Y 50 DIV MUL 2 MUL 1 SUB = intensity.nc
+	gmt makecpt -Chot -T0/10/1
+	gmt subplot begin 4x1 -R0/100/0/50 -Fs15c/5c -JX15c/5c -Scb -Srl -A+gwhite+p0.25p -Bafg10 -X3c -Y3c
+		# 1. Just plot the stripes grid with no other effects
+		gmt subplot set 0 -A"PLAIN GRID and CPT [NO -Q]"
+		gmt basemap
+		echo 50 20 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage stripes.nc
+		# 2. Plot stripes with intensity variations 
+		gmt subplot set 1 -A"SAME PLUS LINEAR INTENSITY FROM LL (-1) TO UR (+1) [No -Q]"
+		gmt basemap
+		echo 50 20 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage stripes.nc -Iintensity.nc
+		# 3. Plot stripes but set z = 5 to be transparent
+		gmt subplot set 2 -A"PLAIN GRID WITH TRANSPARENCY FOR Z = 5 [-Q+z5]"
+		gmt basemap
+		echo 50 20 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage stripes.nc -Q+z5
+		# 4. Plot stripes but set color equal yellow (z = 7-8) to be transparent
+		gmt subplot set 3 -A"PLAIN GRID WITH TRANSPARENCY AT COLOR = yellow (z = 7-8) [-Qyellow]"
+		gmt basemap
+		echo 70 25 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage stripes.nc -Qyellow
+	gmt subplot end
+	# Place the colorbar beneath the subplot */
+	gmt colorbar -DJBC
+gmt end show

test/grdimage/grid_transp.sh

@@ -0,0 +1,35 @@
+#!/usr/bin/env bash
+#
+# Testing various grids and transparency in grdimage
+# Top row is a topo relief grid (image via CPT) and the same grid with fake NaNs to be transparent
+
+gmt begin grid_transp
+	gmt set FONT_TAG 12p PS_MEDIA letter
+	# Make two grids: one a 0|1 multiply grid setting a region to zero and one making NaNs over Africa and 1 elsewhere
+	gmt grdmath -Rd -I01d -rp -fg 1 X -40 SUB ABS 40 LT Y -20 SUB ABS 20 LT MUL SUB = zero.nc
+	gmt grdmath -Rzero.nc 0 0 SDIST 4000 GT 0 NAN = nan.nc
+	gmt subplot begin 4x1 -Rd -Fs12c/6c -JQ12c -Scb -Srl -A+gwhite+p0.25p -Bafg10 -X5c -Y1c
+		# 1. Plot 01d relief grid on top of grid lines (which will be invisible)
+		gmt subplot set 0 -A"TOPOGRAPHY GRID WITH NO NANS"
+		gmt basemap
+		echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage @earth_relief_01d_p
+		# 2. Plot 01d relief grid after setting a R=4000 chunk around (0,0) to NaNs
+		gmt subplot set 1 -A"TOPOGRAPHY GRID WITH NANS AS MISSING DATA"
+		gmt grdmath nan.nc @earth_relief_01d_p MUL = hurt.nc
+		gmt basemap
+		echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage hurt.nc -Cgeo
+		# 3. Plot the same grid but now say NaNs should be transparent
+		gmt subplot set 2 -A"TOPO WITH NANS INDICATING FULL TRANSPARENCY"
+		gmt basemap
+		echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage hurt.nc -Cgeo -Q
+		# 4. Plot 01d relief grid after setting a big chunk to zero and say zero is transparent
+		gmt grdmath zero.nc @earth_relief_01d_p MUL = hurt.nc
+		gmt subplot set 3 -A"TOPO WITH ZEROS INDICATING FULL TRANSPARENCY"
+		gmt basemap
+		echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+		gmt grdimage hurt.nc -Cgeo -Q+z0
+	gmt subplot end
+gmt end show

test/grdimage/image_vartrans.sh

@@ -0,0 +1,31 @@
+#!/usr/bin/env bash
+# For variable transparency TIF we plot pixels as squares with individual transparency
+# Must add -Q to revert to original opaque blended image (color match the transparent
+# colors but the pixels remain opaque).
+
+gmt begin image_vartrans
+	gmt set FONT_TAG 9p PS_MEDIA letter
+	gmt subplot begin 4x1 -Rd -Fs12c/6c -JQ12c -Scb -Srl -A+gwhite+p0.25p -Bafg10 -X5c -Y1c
+	# 1. Create and plot a continuously varying transparency grid (0-1)
+	gmt grdmath -Rd -I01d -r Y 180 DIV 0.5 ADD = transparency.nc
+	gmt subplot set 0 -A"FAKE VARIABLE TRANSPARENCIES"
+	gmt grdcontour transparency.nc -C0.05 -A0.1
+	# 2. Plot original day image
+	gmt subplot set 1 -A"RGB IMAGE WITH NO TRANSPARENCY (OPAQUE)"
+	gmt basemap
+	echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+	gmt grdimage @earth_day_01d
+	# 3. Achieve variable transparency via transparent squares as pixels
+	gmt subplot set 2 -A"RGBA IMAGE WITH VARIABLE TRANSPARENCY AS SQUARES"
+	gmt basemap
+	echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+	# Plot 01d day image grid after adding variable transparency
+	gmt grdmix @earth_day_01d -Atransparency.nc -C -Grgba.tif
+	gmt grdimage rgba.tif
+	# 4. Do variable transparency via image blending instead (opaque)
+	gmt subplot set 3 -A"RGBA IMAGE WITH VARIABLE TRANSPARENCY AS OPAQUE BLEND"
+	gmt basemap
+	echo 0 0 BACKGROUND | gmt text -F+f32p,1+a45
+	gmt grdimage rgba.tif -Q
+	gmt subplot end
+gmt end show

test/grdimage/transp_mix.sh

@@ -0,0 +1,36 @@
+#!/usr/bin/env bash
+#
+# Testing gmt grdmix on adding A to RGB image and plot via grdimage
+
+gmt begin transp_mix ps
+	gmt set MAP_GRID_PEN 2p FONT_TAG 9p PS_MEDIA letter PS_PAGE_ORIENTATION portrait
+	gmt grdmath -Rd -I01d -rp X 360 DIV 0.5 ADD Y 180 DIV 0.5 ADD MUL = transparency.nc
+	gmt subplot begin 3x2 -Rd -Fs9c/4.5c -JQ9c -Scb -Srl -A+gwhite+p0.25p -Bafg30 -X2c
+		# 1. Plot 01d day image on top of grid lines (which will be invisible)
+		gmt subplot set 0 -A"JUST IMAGE (OPAQUE)"
+		gmt basemap
+		gmt grdimage @earth_day_01d_p
+		# 2. Plot the fake transparencies
+		gmt subplot set 1 -A"FAKE TRANSPARENCY (OPAQUE)"
+		gmt basemap
+		gmt grdimage transparency.nc -Chot
+		# 3. Ignore transparencies when making the tif and plot the image on top of grid lines (which will be invisible)
+		gmt grdmix @earth_day_01d_p -C -Grgba.tif -A0.7
+		gmt subplot set 2 -A"IMAGE FIXED 0.7 TRANSPARENCY (TRUE TRANSPARENCY)"
+		gmt basemap
+		gmt grdimage rgba.tif
+		# 4. Mix transparencies into the tif and plot the image on top of grid lines (which will be invisible)
+		gmt grdmix @earth_day_01d_p -Atransparency.nc+o -C -Grgba.tif
+		gmt subplot set 3 -A"BLENDED TRANSPARENT IMAGE (OPAQUE)"
+		gmt basemap
+		gmt grdimage rgba.tif -Q
+		# 5. Plot image using -Q+i to invert opacity to transparency on top of grid lines (which will be invisible)
+		gmt subplot set 4 -A"BLENDED OPACITY IMAGE (OPAQUE)"
+		gmt basemap
+		gmt grdimage rgba.tif -Q
+		# 6. Plot image using defaults on top of grid lines (which will visible through the squares)
+		gmt subplot set 5 -A"DEFAULT IMAGE VIA SQUARES (TRUE TRANSPARENCY)"
+		gmt basemap
+		gmt grdimage rgba.tif
+	gmt subplot end
+gmt end show

test/psscale/interval-panel.sh

@@ -0,0 +1,29 @@
+#!/usr/bin/env bash
+#
+# Testing categorical colorbars with default panels and interval annotations
+# E.g., https://forum.generic-mapping-tools.org/t/colorbar-box-and-annotation/4548
+
+ps=interval-panel.ps
+
+cat << EOF > t.cpt
+# COLOR_MODEL = RGB
+0	247/251/255	10	247/251/255
+10	222/235/247	20	222/235/247
+20	198/219/239	50	198/219/239
+50	158/202/225	100	158/202/225
+100	107/174/214	150	107/174/214
+150	066/146/198	200	066/146/198
+200	033/113/181	300	033/113/181
+300	008/081/156	500	008/081/156
+500	008/048/107	600	008/048/107
+EOF
+# Vertical colorbars
+gmt psscale -Ct.cpt -Dx0c/13c+w12c+v -Li3p -F+p2p -P -K > ${ps}
+gmt psscale -Ct.cpt -Dx5c/13c+w12c+v -LI3p -F+p2p -O -K >> ${ps}
+gmt psscale -Ct.cpt -Dx10c/13c+w12c+v -Li3p -O -K >> ${ps}
+gmt psscale -Ct.cpt -Dx15c/13c+w12c+v -LI3p -O -K >> ${ps}
+# Horizontal colorbars
+gmt psscale -Ct.cpt -Dx0c/1c+w18c+h -Li3p -F+p2p -O -K -X-1c >> ${ps}
+gmt psscale -Ct.cpt -Dx0c/4c+w18c+h -LI3p -F+p2p -O -K >> ${ps}
+gmt psscale -Ct.cpt -Dx0c/7c+w18c+h -Li3p -O -K >> ${ps}
+gmt psscale -Ct.cpt -Dx0c/10c+w18c+h -LI3p -O >> ${ps}