Let grdcut -E extract a vertical slice from a 3-D cube

doc/rst/source/grdcut.rst

@@ -17,6 +17,7 @@ Synopsis
 |-G|\ *outgrid*
 |SYN_OPT-R|
 [ |-D|\ [**+t**] ]
+[ |-E|\ **x**\|\ **y**\ *coord* ]
 [ |-F|\ *polygonfile*\ [**+c**][**+i**] ]
 [ |-J|\ *parameters* ]
 [ |-N|\ [*nodata*] ]
@@ -38,9 +39,10 @@ If in doubt, run :doc:`grdinfo` to check range. Alternatively, define the
 sub-region indirectly via a range check on the node values or via distances from
 a fixed point. Furthermore, you can use |-J| for oblique projections to determine
 the corresponding rectangular |-R| setting that will give a sub-region that fully
-covers the oblique domain.  Finally, use |-F| to specify a polygon and either use
+covers the oblique domain.  You can use |-F| to specify a polygon and either use
 its bounding box for sub-region or set grid nodes inside or outside the polygon
-to NaN. **Note**: If the input grid is actually an image (gray-scale,
+to NaN. Finally, if the input is a 3-D netCDF cube then you can make a vertical
+slice through existing nodes. **Note**: If the input grid is actually an image (gray-scale,
 RGB, or RGBA), then options |-N| and |-Z| are unavailable, while for multi-layer
 Geotiff files only options |-R|, |-S| and |-G| are supported, i.e., you can cut out
 a sub-region only (which we do via *gdal_translate* if you have multiple bands).
@@ -80,6 +82,15 @@ Optional Arguments
     remote gridded data set without implied resolution. Append **+t** to instead receive
     the information as the trailing string "-Rwest/east/south/north -Ixinc/yinc".
 
+.. _-E:
+
+**-E**\ **x**\|\ **y**\ *coord*
+    We extract a vertical slice going along the **x**\ -column *coord* or along the
+    **y**\ -row *coord*, depending on the given directive.
+    **Note**: (1) Input file must be a 3-D netCDF cube, and |-E| can only be used with
+    option |-G|. (2) *coord* must exactly match the coordinates given by th cube. We are
+    not interpolating between nodes and only do a clean slice through existing cube nodes.
+
 .. _-F:
 
 **-F**\ *polygonfile*\ [**+c**][**+i**]
@@ -184,6 +195,11 @@ that is using an oblique projection to display the remote Earth Relief data grid
 
     gmt grdcut @earth_relief -R270/20/305/25+r -JOc280/25.5/22/69/24c -D+t -V
 
+To extract a vertical grid slice at *x = 35* and parallel to the *y-z* plane
+from the 3-D model seis3D.nc, try::
+
+    gmt grdcut seis3D.nc -Gslice_x35.nc -Ex35 -V
+
 Notes
 -----
 

doc/rst/source/module_core_purpose.rst_

@@ -76,7 +76,7 @@
 
 .. |grdconvert_purpose| replace:: Convert between different grid formats
 
-.. |grdcut_purpose| replace:: Extract subregion from a grid or image
+.. |grdcut_purpose| replace:: Extract subregion from a grid or image or a slice from a cube
 
 .. |grdedit_purpose| replace:: Modify header or content of a grid
 

src/gmt_api.c

@@ -6149,6 +6149,8 @@ GMT_LOCAL struct GMT_CUBE * gmtapi_import_cube (struct GMTAPI_CTRL *API, int obj
 					z_min = U_obj->header->z_min;	/* Initialize cube min/max values based on this first layer */
 					z_max = U_obj->header->z_max;
 					here = 0;	/* Initialize offset into k'th layer */
+					UH = gmt_get_U_hidden (U_obj);
+					UH->alloc_mode = GMT_ALLOC_INTERNALLY;
 				}
 				else {	/* Here we update min/max for subsequent layers read */
 					if (G->header->z_min < z_min) z_min = G->header->z_min;
@@ -13421,11 +13423,14 @@ struct GMT_RESOURCE * GMT_Encode_Options (void *V_API, const char *module_name,
 	else if (!strncmp (module, "pscoupe", 7U)) {
 		type = ((opt = GMT_Find_Option (API, 'A', *head)) && strstr (opt->arg, "+c")) ? 'D' : 'X';
 	}
-	/* 1x. Check if grdcut is just getting information and if input is a grid or image */
+	/* 1x. Check if grdcut is just getting information and if input is a grid, image or cube */
 	else if (!strncmp (module, "grdcut", 6U)) {
 		if (GMT_Find_Option (API, 'D', *head))
 			deactivate_output = true;   /* Turn off implicit output since none is in effect, only secondary -D output */
-		type = (GMT_Find_Option (API, 'I', *head)) ? 'I' : 'G'; /* Giving -I means we are reading an image */
+		else if (GMT_Find_Option (API, 'E', *head))
+			type = 'U'; /* Giving -E means we are reading a 3-D cube */
+		else
+			type = 'G'; /* Default is a grid */
 	}
 	/* 1y. Check if this is the gravprisms module, where primary dataset input should be turned off if -C is used */
 	else if (!strncmp (module, "gravprisms", 10U) && (opt = GMT_Find_Option (API, 'C', *head))) {

src/gmt_grdio.c

@@ -3957,3 +3957,74 @@ bool gmt_grd_domains_match (struct GMT_CTRL *GMT, struct GMT_GRID *A, struct GMT
 	}
 	return (true);
 }
+
+struct GMT_GRID * gmt_vertical_cube_cut (struct GMT_CTRL *GMT, struct GMT_CUBE *C, unsigned int dim, double coord) {
+	/* Special case of slicing a cube vertically and along the cube's node structure.  For
+	 * oblique cuts and resampling and for grids with variable spacing/time in the third
+	 * dimension, see grdinterpolate instead.
+	 * Here, dim is either GMT_X or GMT_Y which informs us whether coord is an x-
+	 * or y-coordinate and that defines the vertical plan to be at that constant
+	 * coordinate and parallel to the z-axis and the other axis (y or x).
+	 */
+	uint64_t row, col, xrow, xcol, layer, ijg, ijc;
+	double pos = 0.0;
+	struct GMT_GRID *G = NULL;
+	struct GMT_GRID_HIDDEN *GH = NULL;
+
+	if (gmtlib_var_inc (C->z, C->header->n_bands)) {	/* Check if equidistant in z direction */
+		GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cube has non-equidistant spacing in the third dimension (consider grdinterpolate instead)\n");
+		GMT->parent->error = GMT_RUNTIME_ERROR;
+		return (NULL);
+	}
+	G = gmt_create_grid (GMT);	/* Create empty grid structure */
+	/* The number of columns in the output grid depends which of the two planes we selected */
+	G->header->n_columns  = (dim == GMT_X) ? C->header->n_rows     : C->header->n_columns;
+	G->header->wesn[XLO]  = (dim == GMT_X) ? C->header->wesn[YLO]  : C->header->wesn[XLO];
+	G->header->wesn[XHI]  = (dim == GMT_X) ? C->header->wesn[YHI]  : C->header->wesn[XHI];
+	G->header->inc[GMT_X] = (dim == GMT_X) ? C->header->inc[GMT_Y] : C->header->inc[GMT_X];
+	/* The number of rows in the output grid is always the third dimension in the cube and independent of dim */
+	G->header->n_rows     = C->header->n_bands;	
+	G->header->wesn[YLO]  = C->z_range[0];
+	G->header->wesn[YHI]  = C->z_range[1];
+	G->header->inc[GMT_Y] = C->z_inc;
+	gmt_set_grddim (GMT, G->header);	/* Determine dimensions, mx, pad etc */
+	if (dim == GMT_X) {	/* Received an x-coordinate, find corresponding column */
+		xcol = gmt_M_grd_x_to_col (GMT, coord, C->header);
+		pos = gmt_M_grd_col_to_x (GMT, xcol, C->header);
+	}
+	else {	/* Received an y-coordinate, find corresponding row */
+		xrow = gmt_M_grd_y_to_row (GMT, coord, C->header);
+		pos = gmt_M_grd_row_to_y (GMT, xrow, C->header);
+	}
+	if (!doubleAlmostEqualZero (coord, pos)) {	/* Not aligned with cube nodes */
+		static char *axis = "xy";
+		GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -E%c: Your %c-coordinate is not aligned with the cube %c-nodes (%lg vs %lg).\n", axis[dim], axis[dim], axis[dim], coord, pos);
+		GMT->parent->error = GMT_RUNTIME_ERROR;
+		gmt_free_grid (GMT, &G, true);
+		return (NULL);
+	}
+
+	/* Checks passed, time to do the work */
+
+	G->data = gmt_M_memory_aligned (GMT, NULL, G->header->size, gmt_grdfloat);	/* Allocate grid and padding */
+	GH = gmt_get_G_hidden (G);
+	GH->alloc_mode = GMT_ALLOC_INTERNALLY;
+
+	/* Loop over the output grids rows and columns and match to nodes in the cube */
+	for (layer = 0; layer < G->header->n_rows; layer++) {	/* This is the loop over the rows in the output grid */
+		ijg = gmt_M_ijp (G->header, G->header->n_rows-1-layer, 0);	/* TL node in output grid */
+		if (dim == GMT_X) {	/* Must loop over the cube's y-dimension */
+			ijc = gmt_M_ijp (C->header, 0, xcol) + layer * C->header->size;	/* Corresponding point in the cube */
+			/* col loops over the columns in the 2-D grid which are rows in the cube */
+			for (col = 0; col < G->header->n_columns; col++, ijg++, ijc += G->header->mx)
+				G->data[ijg] = C->data[ijc];
+		}
+		else {	/* Must loop over the cube's x-dimension */
+			ijc = gmt_M_ijp (C->header, xrow, 0) + layer * C->header->size;/* Corresponding point in the cube */
+			/* col here is the column in the 2-D grid which are also cols in the cube */
+			for (col = 0; col < C->header->n_columns; col++, ijg++, ijc++)
+				G->data[ijg] = C->data[ijc];
+		}
+	}
+	return (G);	/* Return the fully allocated local grid */
+}

src/gmt_prototypes.h

@@ -218,6 +218,7 @@ EXTERN_MSC bool gmt_file_is_cache (struct GMTAPI_CTRL *API, const char *file);
 
 /* gmt_grdio.c: */
 
+EXTERN_MSC struct GMT_GRID * gmt_vertical_cube_cut (struct GMT_CTRL *GMT, struct GMT_CUBE *C, unsigned int dim, double coord);
 EXTERN_MSC bool gmt_grd_domains_match (struct GMT_CTRL *GMT, struct GMT_GRID *A, struct GMT_GRID *B, char *comment);
 EXTERN_MSC unsigned int gmt_grid_perimeter (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h, double **x, double **y);
 EXTERN_MSC void gmt_change_grid_history (struct GMTAPI_CTRL *API, unsigned int mode, struct GMT_GRID_HEADER *h, char *command);

src/grdclip.c

@@ -33,11 +33,11 @@
 
 #define THIS_MODULE_CLASSIC_NAME	"grdclip"
 #define THIS_MODULE_MODERN_NAME		"grdclip"
-#define THIS_MODULE_LIB			"core"
-#define THIS_MODULE_PURPOSE		"Clip the range of grid values"
-#define THIS_MODULE_KEYS		"<G{,GG}"
-#define THIS_MODULE_NEEDS		""
-#define THIS_MODULE_OPTIONS 		"-RV"
+#define THIS_MODULE_LIB		"core"
+#define THIS_MODULE_PURPOSE	"Clip the range of grid values"
+#define THIS_MODULE_KEYS	"<G{,GG}"
+#define THIS_MODULE_NEEDS	""
+#define THIS_MODULE_OPTIONS 	"-RV"
 
 enum Grdclip_cases {
 	GRDCLIP_BELOW	= 1,

src/grdcut.c

@@ -32,8 +32,8 @@
 #define THIS_MODULE_CLASSIC_NAME	"grdcut"
 #define THIS_MODULE_MODERN_NAME	"grdcut"
 #define THIS_MODULE_LIB		"core"
-#define THIS_MODULE_PURPOSE	"Extract subregion from a grid or image"
-#define THIS_MODULE_KEYS	"<G{,FD(=,>DD,G?}"
+#define THIS_MODULE_PURPOSE	"Extract subregion from a grid or image or a slice from a cube"
+#define THIS_MODULE_KEYS	"<?{,FD(=,>DD,G?}"
 #define THIS_MODULE_NEEDS	""
 #define THIS_MODULE_OPTIONS "-JRVf"
 
@@ -56,6 +56,11 @@ struct GRDCUT_CTRL {
 		bool text;
 		bool quit;
 	} D;
+	struct GRDCUT_E {	/* -Dx|y<value> */
+		bool active;
+		unsigned dim;	/* GMT_X or GMT_Y */
+		double coord;	/* THe x- or y-value we wish to cut a cube vertically */
+	} E;
 	struct GRDCUT_F {	/* -Fpolfile[+c][+i] */
 		bool active;
 		bool crop;
@@ -114,7 +119,7 @@ static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *C) {	/* Dealloc
 static int usage (struct GMTAPI_CTRL *API, int level) {
 	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
 	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
-	GMT_Usage (API, 0, "usage: %s %s -G%s %s [-D[+t]] [-F<polygontable>[+c][+i]] [%s] [-N[<nodata>]] [-S<lon>/<lat>/<radius>[+n]] [%s] [-Z[<min>/<max>][+n|N|r]] [%s] [%s]\n",
+	GMT_Usage (API, 0, "usage: %s %s -G%s %s [-D[+t]] [-Ex|y<coord>] [-F<polygontable>[+c][+i]] [%s] [-N[<nodata>]] [-S<lon>/<lat>/<radius>[+n]] [%s] [-Z[<min>/<max>][+n|N|r]] [%s] [%s]\n",
 		name, GMT_INGRID, GMT_OUTGRID, GMT_Rgeo_OPT, GMT_J_OPT, GMT_V_OPT, GMT_f_OPT, GMT_PAR_OPT);
 
 	if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
@@ -130,6 +135,8 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
 	GMT_Usage (API, 1, "\n-D[+t]");
 	GMT_Usage (API, -2, "Dry-run mode. No grid is written but its domain and increment will be "
 		"written to standard output in w e s n dx dy numerical format.  Append +t to instead receive text strings -Rw/e/s/n -Idx/dy.");
+	GMT_Usage (API, 1, "\n-Ex|y<coord>");
+	GMT_Usage (API, -2, "Cut a vertical grid from an input cube along x = <coord> or y == <coord>.\n");
 	GMT_Usage (API, 1, "\n-F<polygontable>[+c][+i]");
 	GMT_Usage (API, -2, "Specify a multi-segment closed polygon table that describes the grid subset "
 		"to be extracted (nodes between grid boundary and polygons will be set to NaN).");
@@ -164,7 +171,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT
 	 */
 
 	bool F_or_R_or_J, do_file_check = true;
-	unsigned int n_errors = 0, k;
+	unsigned int n_errors = 0, k, n_files = 0;
 	char za[GMT_LEN64] = {""}, zb[GMT_LEN64] = {""}, zc[GMT_LEN64] = {""}, *c = NULL;
 	struct GMT_OPTION *opt = NULL;
 	struct GMTAPI_CTRL *API = GMT->parent;
@@ -179,6 +186,7 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->In.active);
 				n_errors += gmt_get_required_string (GMT, opt->arg, opt->option, 0, &(Ctrl->In.file));
 				if (do_file_check && GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
+				n_files++;
 			break;
 
 			/* Processes program-specific parameters */
@@ -191,6 +199,22 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT
 					gmt_M_str_free (opt->arg);	/* Free internal marker */
 				}
 				break;
+			case 'E':
+				n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
+				switch (opt->arg[0]) {
+					case 'x':	Ctrl->E.dim = GMT_X;	break;
+					case 'y':	Ctrl->E.dim = GMT_Y;	break;
+					default:
+						GMT_Report (API, GMT_MSG_ERROR, "Option -E: Must select directive x or y\n");
+						n_errors++;
+						break;
+				}
+				if (opt->arg[1]) {
+					n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, Ctrl->E.dim),
+					                                     gmt_scanf_arg (GMT, &opt->arg[1], gmt_M_type (GMT, GMT_IN, Ctrl->E.dim), false,
+					                                     &Ctrl->E.coord), &opt->arg[1]);
+				}
+				break;
 			case 'F':
 				n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
 				if (opt->arg[0]) {
@@ -286,14 +310,18 @@ static int parse (struct GMT_CTRL *GMT, struct GRDCUT_CTRL *Ctrl, struct GMT_OPT
 	}
 
 	n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && Ctrl->G.file, "Option -D: Cannot specify -G since no grid will be returned\n");
-	//n_errors += gmt_M_check_condition (GMT, Ctrl->D.active && !GMT->common.J.active, "Option -D: Requires -R and -J\n");
 	n_errors += gmt_M_check_condition (GMT, GMT->common.R.active[RSET] && Ctrl->F.crop, "Option -F: Modifier +c cannot be used with -R\n");
 	F_or_R_or_J = GMT->common.R.active[RSET] || Ctrl->F.active || GMT->common.J.active;
-	n_errors += gmt_M_check_condition (GMT, (F_or_R_or_J + Ctrl->S.active + Ctrl->Z.active) != 1,
+	n_errors += gmt_M_check_condition (GMT, (F_or_R_or_J + Ctrl->S.active + Ctrl->E.active+ Ctrl->Z.active) != 1,
 	                                   "Must specify only one of the -F, -R, -S or the -Z options\n");
 	n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file && !Ctrl->D.active, "Option -G: Must specify output grid file\n");
 	n_errors += gmt_M_check_condition (GMT, !Ctrl->In.active, "Must specify one input grid file\n");
-	if (n_errors == 0) {
+	n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && !Ctrl->G.active, "Option -E: Must specify output grid file\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && n_files != 1, "Option -E: Must supply an input cube\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && !gmt_nc_is_cube (API, Ctrl->In.file), "Option -E: Must supply an input cube, not grid\n");
+	n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && (F_or_R_or_J + Ctrl->N.active + Ctrl->S.active + Ctrl->Z.active), "Option -E: Can only be used with -G\n");
+
+	if (n_errors == 0 && !Ctrl->E.active) {
 		if (!Ctrl->D.quit) {
 			int ftype = gmt_raster_type(GMT, Ctrl->In.file, true);
 			if (ftype == GMT_IS_IMAGE)	/* Must read file as an image */
@@ -416,7 +444,6 @@ GMT_LOCAL unsigned int grdcut_node_is_outside (struct GMT_CTRL *GMT, struct GMT_
 	return ((inside) ? 0 : 1);	/* 1 if outside */
 }
 
-
 #define bailout(code) {gmt_M_free_options (mode); return (code);}
 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
 
@@ -461,7 +488,22 @@ EXTERN_MSC int GMT_grdcut (void *V_API, int mode, void *args) {
 	if (Ctrl->D.quit)	/* Already reported information deep inside gmt_init_module */
 		Return (GMT_NOERROR);
 
-	gmt_grd_set_datapadding (GMT, true);	/* Turn on gridpadding when reading a subset */
+	gmt_grd_set_datapadding (GMT, true);	/* Turn on grid padding when reading a subset */
+
+	if (Ctrl->E.active) {	/* Extract a vertical slice grid aligned with x or y axis from a cube */
+		struct GMT_CUBE *C = NULL;
+		if ((C = GMT_Read_Data (API, GMT_IS_CUBE, GMT_IS_FILE, GMT_IS_VOLUME, GMT_CONTAINER_AND_DATA, NULL, Ctrl->In.file, NULL)) == NULL)
+			Return (GMT_DATA_READ_ERROR);
+		if ((G = gmt_vertical_cube_cut (GMT, C, Ctrl->E.dim, Ctrl->E.coord)) == NULL)
+			Return (API->error);
+		if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, G))
+			Return (API->error);
+		if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, G) != GMT_NOERROR) {
+			Return (API->error);
+		}
+		gmt_free_grid (GMT, &G, true);
+		Return (GMT_NOERROR);
+	}
 
 	if (!Ctrl->Z.active) {	/* All of -F, -R, -S selections first needs the header */
 		if (Ctrl->In.type == GMT_IS_IMAGE) {

test/baseline/grdcut.dvc

@@ -1,5 +1,5 @@
 outs:
-- md5: 9b894ca37d8aebb68012a7076f065cae.dir
-  size: 566243
-  nfiles: 5
+- md5: d0bfce2251b06e3137263b75eaacf6c0.dir
+  size: 736923
+  nfiles: 6
   path: grdcut

test/grdcut/vert_cube_cuts.sh

@@ -0,0 +1,25 @@
+#!/usr/bin/env bash
+# Testing vertical cuts of a cube via grdcut
+#
+gmt begin vert_cube_cuts
+	# Create the 3-D cube from w = 5 * exp (-0.5*(r/z)^4) + 8 - z
+	gmt math -T1/8/1 -o0 --FORMAT_FLOAT_OUT=%01.0f T = z.txt
+	while read z; do
+		gmt grdmath -R-10/10/-10/10 -I1 0 0 CDIST $z 0.25 MUL HYPOT $z DIV 4 POW -0.5 MUL EXP 5 MUL 8 $z SUB ADD = tmp_$z.grd
+	done < z.txt
+	gmt grdinterpolate -Zz.txt tmp_?.grd -Gfake_xyz_cube.nc
+	# Cut slices at x = 0 and y = 6
+	gmt grdcut fake_xyz_cube.nc -Ex0 -Gy.grd
+	gmt grdcut fake_xyz_cube.nc -Ey6 -Gx.grd
+	# Compute the yx and xz vertical grids directly using -R with x|y and z
+	gmt grdmath -R-10/10/1/8 -I1 X Y 0.25 MUL HYPOT Y DIV 4 POW -0.5 MUL EXP 5 MUL 8 Y SUB ADD = orig_y.grd
+	gmt grdmath -R-10/10/1/8 -I1 X 6 HYPOT Y 0.25 MUL HYPOT Y DIV 4 POW -0.5 MUL EXP 5 MUL 8 Y SUB ADD = orig_x.grd
+	gmt makecpt -T0/12
+	gmt subplot begin 2x2 -Fs8c -Sct+tc -Srl -Blrbt -A1+gwhite -R-10/10/1/8 -JX10c/-5c -M0.5c -T"Cuts of @[5 e^{\left [-\frac{1}{2}\left( \frac{r}{z}\right)^4\right]} + 8 - z@["
+		gmt grdimage orig_x.grd -c -B+t"y-z cut x = 0"
+		gmt grdimage orig_y.grd -c -B+t"x-z cut y = 6"
+		gmt grdimage x.grd -c
+		gmt grdimage y.grd -c
+	gmt subplot end
+	gmt colorbar -DJBC -B
+gmt end show