Enable writing of 3-D netCDF data cubes from greenspline and grdinterpolate

src/gmt_api.c

@@ -1997,6 +1997,11 @@ GMT_LOCAL int gmtapi_init_grdheader (struct GMT_CTRL *GMT, unsigned int directio
 	gmtlib_grd_get_units (GMT, header);
 	gmt_BC_init (GMT, header);	/* Initialize grid interpolation and boundary condition parameters */
 	HH->grdtype = gmtlib_get_grdtype (GMT, direction, header);	/* Set grid type (i.e. periodicity for global grids) */
+#ifdef DOUBLE_PRECISION_GRID
+	header->type = GMT_GRID_IS_ND;
+#else
+	header->type = GMT_GRID_IS_NF;
+#endif
 	return (GMT_NOERROR);
 }
 
@@ -5372,7 +5377,7 @@ GMT_LOCAL struct GMT_CUBE * gmtapi_import_cube (struct GMTAPI_CTRL *API, int obj
 	unsigned int both_set = (GMT_CONTAINER_ONLY | GMT_DATA_ONLY);
 	unsigned int method, start_over_method = 0;
 	double dx, dy, d;
-	double *level = NULL, z_min, z_max;
+	double *level = NULL, z_min, z_max, w_range[2] = {0.0, 0.0};
 	p_func_uint64_t GMT_2D_to_index = NULL;
 	struct GMT_CUBE *U_obj = NULL, *U_orig = NULL;
 	struct GMT_GRID *G = NULL;
@@ -5435,7 +5440,7 @@ GMT_LOCAL struct GMT_CUBE * gmtapi_import_cube (struct GMTAPI_CTRL *API, int obj
 						strcpy (cube_layer, &nc_z_named[1]);	/* Place variable name in cube_layer string */
 						nc_z_named[0] = '\0';	/* Chop off layer name for now */
 					}
-					if (gmt_examine_nc_cube (GMT, the_file, &n_layers, &level)) {	/* Learn the basics about the cube */
+					if (gmt_nc_read_cube_info (GMT, the_file, w_range, &n_layers, &level)) {	/* Learn the basics about the cube */
 						GMT_Report (API, GMT_MSG_ERROR, "gmtapi_import_cube: Unable to examine cube %s.\n", the_file);
 						return_null (API, GMT_RUNTIME_ERROR);
 					}
@@ -5880,14 +5885,13 @@ GMT_LOCAL int gmtapi_export_cube (struct GMTAPI_CTRL *API, int object_ID, unsign
 	bool done = true;
 	unsigned int method;
 	uint64_t row, col, i0, i1, j0, j1, k0, k1, ij, ijp, ij_orig;
-	uint64_t k, n_layers_used, here = 0, save_n_bands;
+	uint64_t k, here = 0;
 	size_t size;
 	double dx, dy;
 	p_func_uint64_t GMT_2D_to_index = NULL;
 	GMT_putfunction api_put_val = NULL;
 	struct GMTAPI_DATA_OBJECT *S_obj = NULL;
 	struct GMT_CUBE *U_copy = NULL;
-	struct GMT_GRID *G = NULL;
 	struct GMT_MATRIX *M_obj = NULL;
 	struct GMT_MATRIX_HIDDEN *MH = NULL;
 	struct GMT_CUBE_HIDDEN *UH = gmt_get_U_hidden (U_obj), *UH2 = NULL;
@@ -5911,52 +5915,18 @@ GMT_LOCAL int gmtapi_export_cube (struct GMTAPI_CTRL *API, int object_ID, unsign
 	}
 	if (mode & GMT_GRID_IS_GEO) gmt_set_geographic (GMT, GMT_OUT);	/* From API to tell cube is geographic */
 	gmtlib_grd_set_units (GMT, U_obj->header);	/* Ensure unit strings are set, regardless of destination */
+
 	method = gmtapi_set_method (S_obj);	/* Get the actual method to use since may be MATRIX or VECTOR masquerading as GRID */
 	switch (method) {
-		case GMT_IS_FILE:	/* Name of a cube file on disk */
+		case GMT_IS_FILE:	/* Name of a cube file to write to disk */
 			if (mode & GMT_CONTAINER_ONLY) {	/* Update header structure only */
 				GMT_Report (API, GMT_MSG_INFORMATION, "Updating cube header for file %s not implemented\n", S_obj->filename);
 				return (gmtlib_report_error (API, GMT_RUNTIME_ERROR));
 			}
 			else {
 				GMT_Report (API, GMT_MSG_INFORMATION, "Writing cube to file %s\n", S_obj->filename);
-				/* Determine which layers we want to write */
-				k0 = 0;	k1 = U_obj->header->n_bands - 1;
-				if (S_obj->region && S_obj->wesn[ZHI] > S_obj->wesn[ZLO]) {	/* Want to write a subset of layers */
-					if (gmt_get_active_layers (GMT, U_obj, &(S_obj->wesn[ZLO]), &k0, &k1) == 0) {
-						gmtlib_report_error (API, GMT_RUNTIME_ERROR);
-					}
-				}
-				n_layers_used = k1 - k0 + 1;	/* Total number of layers actually to be read */
-				if (n_layers_used == 0) {
-					GMT_Report (API, GMT_MSG_ERROR, "gmtapi_export_cube: No layers selected from GMT_IS_CUBE.\n");
-					return (gmtlib_report_error (API, GMT_DIM_TOO_SMALL));
-				}
-
-				/* !! Remember that the gmt_nc.c codes will unpad and mess up C->data */
-				G = gmt_get_grid (GMT);	/* Need a dummy grid for writing */
-				save_n_bands = U_obj->header->n_bands;	/* Remember how many bands */
-				G->header = U_obj->header;	/* Use this pointer for now */
-				G->header->n_bands = 1;	/* Grids only have one band */
-				here = k0 * U_obj->header->size;	/* Start position in the cube for layer k0 */
-				/* Write the layers via a grid */
-				for (k = k0; k <= k1; k++) {	/* For all selected output levels */
-					if (n_layers_used > 1)	/* Create the k'th layer file */
-						sprintf (file, S_obj->filename, U_obj->z[k]);
-					else	/* Just this one layer grid */
-						sprintf (file, "%s", S_obj->filename);
-					G->data = &U_obj->data[here];	/* Point to start of this layer */
-					GMT_Report (API, GMT_MSG_DEBUG, "gmtapi_export_cube: Layer %" PRIu64 ", offset = %" PRIu64 ".\n", k, here);
-					if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, S_obj->wesn, file, G) != GMT_NOERROR) {
-						return (API->error);
-					}
-					here += U_obj->header->size;	/* Move to next level */
-				}
-				U_obj->header->n_bands = save_n_bands;	/* Restore number of layers */
-				/* Wipe and free the temporary grid structure */
-				G->data = NULL;
-				G->header = NULL;
-				gmt_free_grid (GMT, &G, true);
+				if (gmt_nc_write_cube (GMT, U_obj, S_obj->wesn, S_obj->filename) != GMT_NOERROR)
+					return (gmtlib_report_error (API, API->error));
 				done = true;
 			}
 			break;
@@ -12453,7 +12423,13 @@ struct GMT_RESOURCE * GMT_Encode_Options (void *V_API, const char *module_name,
 	else if (!strncmp (module, "greenspline", 11U) && (opt = GMT_Find_Option (API, 'R', *head))) {
 		/* Found the -R"domain" option; determine the dimensionality of the output */
 		unsigned dim = gmtapi_determine_dimension (API, opt->arg);
-		type = (dim == 2) ? 'G' : 'D';
+		switch (dim) {	/* Determine if output is D, G, or U */
+			case 1: type = 'D'; break;	/* 1-D is a data table */
+			case 2: type = 'G'; break;	/* 2-D is always a grid */
+			default:	/* 3-D, but can be dataset or cube */
+				type = ((opt = GMT_Find_Option (API, 'G', *head))) ? 'U' : 'D';
+				break;
+		}
 	}
 	/* 1g. Check if this is the triangulate module, where primary dataset output should be turned off if -G given without -M,N,Q,S */
 	else if (!strncmp (module, "triangulate", 11U) && (opt = GMT_Find_Option (API, 'G', *head))) {

src/gmt_grdio.c

@@ -2692,6 +2692,11 @@ struct GMT_GRID *gmt_create_grid (struct GMT_CTRL *GMT) {
 	G->header = gmt_get_header (GMT);
 	gmt_grd_init (GMT, G->header, NULL, false); /* Set default values */
 	GMT_Set_Index (GMT->parent, G->header, GMT_GRID_LAYOUT);
+#ifdef DOUBLE_PRECISION_GRID
+	G->header->type = GMT_GRID_IS_ND;
+#else
+	G->header->type = GMT_GRID_IS_NF;
+#endif
 	GH->alloc_mode = GMT_ALLOC_INTERNALLY;		/* Memory can be freed by GMT. */
 	GH->alloc_level = GMT->hidden.func_level;	/* Must be freed at this level. */
 	GH->id = GMT->parent->unique_var_ID++;		/* Give unique identifier */
@@ -3330,6 +3335,11 @@ struct GMT_CUBE *gmtlib_create_cube (struct GMT_CTRL *GMT) {
 	GU = gmt_get_U_hidden (C);
 	C->header = gmt_get_header (GMT);
 	gmt_grd_init (GMT, C->header, NULL, false); /* Set default values */
+#ifdef DOUBLE_PRECISION_GRID
+	C->header->type = GMT_GRID_IS_ND;
+#else
+	C->header->type = GMT_GRID_IS_NF;
+#endif
 	GMT_Set_Index (GMT->parent, C->header, GMT_GRID_LAYOUT);
 	GU->alloc_mode = GMT_ALLOC_INTERNALLY;		/* Memory can be freed by GMT. */
 	GU->alloc_level = GMT->hidden.func_level;	/* Must be freed at this level. */

src/gmt_hidden.h

@@ -156,6 +156,7 @@ struct GMT_GRID_HEADER_HIDDEN {
 	int z_id;                        /* NetCDF: id of z field */
 	int ncid;                        /* NetCDF: file ID */
 	int xy_dim[2];                   /* NetCDF: dimension order of x and y; normally {1, 0} */
+	int xyz_id[3];                   /* NetCDF: id of x, y, and z (if cube) field */
 	size_t t_index[3];               /* NetCDF: index of higher coordinates */
 	size_t data_offset;              /* NetCDF: distance from the beginning of the in-memory grid */
 	size_t n_alloc;                  /* Bytes allocated for this grid */