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*] ]
@@ -80,6 +81,13 @@ 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**: Input file must be a 3-D netCDF cube, and |-E| can only be used with option |-G|.
+
 .. _-F:
 
 **-F**\ *polygonfile*\ [**+c**][**+i**]

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

@@ -13421,11 +13421,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,47 @@ 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) {
+	uint64_t row, col, xrow, xcol, layer, ijg, ijc;
+	struct GMT_GRID *G = NULL;
+
+	if (gmtlib_var_inc (C->z, C->header->n_bands)) {
+		GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cube has non-equidistant spacing in the third dimension\n");
+		GMT->parent->error = GMT_RUNTIME_ERROR;
+		return (NULL);
+	}
+	G = gmt_create_grid (GMT);
+	G->header->n_columns = C->header->n_bands;	/* Vertical dimension */
+	G->header->n_rows = (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];
+	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);
+	if (dim == GMT_X)
+		xcol = gmt_M_grd_x_to_col (GMT, coord, C->header);
+	else
+		xrow = gmt_M_grd_y_to_row (GMT, coord, C->header);
+	G->data = gmt_M_memory_aligned (GMT, NULL, G->header->size, gmt_grdfloat);
+
+	for (layer = 0; layer < C->header->n_bands; layer++) {
+		ijg = gmt_M_ijp (G->header, layer, 0);
+		if (dim == GMT_X) {
+			ijc = gmt_M_ijp (C->header, 0, xcol) + layer * C->header->size;
+			/* row here is basically the column in the 2-D grid */
+			for (row = 0; row < C->header->n_rows; row++, ijg++)
+				G->data[ijg] = C->data[ijc];
+		}
+		else {
+			ijc = gmt_M_ijp (C->header, xrow, 0) + layer * C->header->size;
+			/* col here is the column in the 2-D grid */
+			ijc = gmt_M_ijp (C->header, xrow, 0) + layer * C->header->size;
+			for (col = 0; col < C->header->n_columns; col++, ijg++)
+				G->data[ijg] = C->data[ijc];
+		}
+	}
+	return (G);
+}

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_Destroy_Data (API, &C);
+		Return (GMT_NOERROR);
+	}
 
 	if (!Ctrl->Z.active) {	/* All of -F, -R, -S selections first needs the header */
 		if (Ctrl->In.type == GMT_IS_IMAGE) {