WIP Implement support for 3-D grids

src/gmt.h

@@ -103,6 +103,7 @@ EXTERN_MSC void * GMT_Get_Vector       (void *API, struct GMT_VECTOR *V, unsigne
 EXTERN_MSC void * GMT_Get_Matrix       (void *API, struct GMT_MATRIX *M);
 EXTERN_MSC int GMT_Put_Vector          (void *API, struct GMT_VECTOR *V, unsigned int col, unsigned int type, void *vector);
 EXTERN_MSC int GMT_Put_Matrix          (void *API, struct GMT_MATRIX *M, unsigned int type, int pad, void *matrix);
+EXTERN_MSC int GMT_Put_Layers          (void *API, struct GMT_GRID_HEADER *header, unsigned int n_layers, double *layer);
 /* These 2 functions are new in 6.0 and are being considered beta */
 EXTERN_MSC int GMT_Put_Strings         (void *API, unsigned int family, void *object, char **array);
 EXTERN_MSC char ** GMT_Get_Strings     (void *API, unsigned int family, void *object);
@@ -113,6 +114,7 @@ EXTERN_MSC uint64_t GMT_Get_Index      (void *API, struct GMT_GRID_HEADER *heade
 EXTERN_MSC double * GMT_Get_Coord      (void *API, unsigned int family, unsigned int dim, void *container);
 EXTERN_MSC int GMT_Set_Index	       (void *API, struct GMT_GRID_HEADER *header, char *code);	/* Experimental */
 EXTERN_MSC uint64_t GMT_Get_Pixel      (void *API, struct GMT_GRID_HEADER *header, int row, int col, int layer);	/* Experimental */
+EXTERN_MSC uint64_t GMT_Get_Index3     (void *API, struct GMT_GRID_HEADER *header, int row, int col, int layer);
 
 /* 11 functions to show and inquire about GMT common options, GMT default settings, object metadata, convert strings to doubles, and message and report printing */
 

src/gmt_api.c

@@ -10018,6 +10018,28 @@ int GMT_Get_Info_ (unsigned int *family, void *container, unsigned int *geometry
 }
 #endif
 
+int GMT_Put_Layers (void *V_API, struct GMT_GRID_HEADER *header, unsigned int n_layers, double *layer) {
+	/*Set the grid header's internal layer information */
+	struct GMTAPI_CTRL *API = NULL;
+
+	if (V_API == NULL) return_error (V_API, GMT_NOT_A_SESSION);
+	if (header == NULL) return_error (API, GMT_PTR_IS_NULL);	/* Error if data is NULL */
+	if (layer == NULL) return_error (API, GMT_PTR_IS_NULL);	/* Error if data is NULL */
+	if (n_layers == 0) return_error (API, GMT_VALUE_NOT_SET);	/* Error if no layers */
+
+	API = gmtapi_get_api_ptr (V_API);
+	API->error = GMT_NOERROR;
+	header->n_layers = n_layers;
+	if (header->layer) gmt_M_free (API->GMT, header->layer);	/* Free whatever was there */
+	gmt_M_memcpy (header->layer, layer, n_layers, double);		/* Copy over the layer values */
+	return (API->error);
+}
+
+#ifdef FORTRAN_API
+int GMT_Put_Layers_ (struct GMT_GRID_HEADER *header, unsigned int *n_layers, double *layer) {
+	return (GMT_Put_Layers (GMT_FORTRAN, header, *n_layers, layer));
+#endif
+
 /*! Convenience function to get grid or image node */
 uint64_t GMT_Get_Index (void *V_API, struct GMT_GRID_HEADER *header, int row, int col) {
 	/* V_API not used but all API functions take V_API so no exceptions! */
@@ -10047,6 +10069,21 @@ uint64_t GMT_Get_Pixel_ (void *h, int *row, int *col, int *layer) {
 }
 #endif
 
+/*! Convenience function to get 3-D grid node */
+uint64_t GMT_Get_Index3 (void *V_API, struct GMT_GRID_HEADER *header, int row, int col, int layer) {
+	/* V_API not used but all API functions take V_API so no exceptions! */
+	gmt_M_unused(V_API);
+	return (GMTAPI_index_function (header, row, col, 0) * layer);
+}
+
+#ifdef FORTRAN_API
+uint64_t GMT_Get_Index3_ (void *h, int *row, int *col, int *layer) {
+	/* Fortran version: We pass the global GMT_FORTRAN structure */
+	return (GMT_Get_Index3 (GMT_FORTRAN, h, *row, *col, *layer));
+}
+#endif
+
+
 /*! Specify image memory layout */
 int GMT_Set_Index (void *V_API, struct GMT_GRID_HEADER *header, char *code) {
 	struct GMTAPI_CTRL *API = NULL;

src/gmt_common.h

@@ -108,9 +108,10 @@ struct GMT_COMMON {
 		uint32_t registration;	/* Registration mode of a grid given via -r or -Rgrid */
 		int row_order;	/* Order of rows in NetCDF output: 0 (not set) or k_nc_start_north or k_nc_start_south */
 		unsigned int mode;	/* For modern mode only: 0 = get exact region from data, 1 = rounded region from data */
+		unsigned int dimension;	/* 2 for 2-D grid regions, 3 for 3-D regions */
 		double wesn[6];		/* Boundaries of west, east, south, north, low-z and hi-z */
 		double wesn_orig[4];	/* Original Boundaries of west, east, south, north (oblique projection may reset wesn above) */
-		double inc[2];	/* For grid increments set via -Idx/dy or implicitly via -Ggrid */
+		double inc[3];	/* For grid increments set via -Idx/dy[/dz] or implicitly via -Ggrid */
 		char string[GMT_LEN256];
 	} R;
 	struct U {	/* -U */

src/gmt_grd.h

@@ -151,7 +151,12 @@ enum gmt_enum_wesnids {
 #define gmt_M_ij(h,row,col) ((uint64_t)(((int64_t)(row))*((int64_t)h->mx)+(int64_t)(col)))
 /*! IJPGI macro using h and the pad info that works for either grids (n_bands = 1) or images (n_bands = 1,3,4) */
 #define gmt_M_ijpgi(h,row,col) ((uint64_t)(((int64_t)(row)+(int64_t)h->pad[YHI])*((int64_t)h->mx*(int64_t)h->n_bands)+((int64_t)(col)+(int64_t)h->pad[XLO])*(int64_t)h->n_bands))
-
+/*! IJKP macro using h and the pad info */
+#define gmt_M_ijkp(h,row,col,layer) (((uint64_t)layer) * gmt_M_ijp(h,row,col))
+/*! IJ0 macro using h but ignores the pad info */
+#define gmt_M_ijk0(h,row,col,layer) (((uint64_t)layer) * gmt_M_ij0(h,row,col))
+/*! IJ macro using h but treats the entire grid with pad as no-pad grid, i.e. using mx as width */
+#define gmt_M_ijk(h,row,col,layer) (((uint64_t)layer) * gmt_M_ij(h,row,col))
 /*! Obtain row and col from index */
 #define gmt_M_col(h,ij) (((ij) % h->mx) - h->pad[XLO])
 #define gmt_M_row(h,ij) (((ij) / h->mx) - h->pad[YHI])

src/gmt_grdio.c

@@ -100,7 +100,7 @@ GMT_LOCAL inline struct GMT_IMAGE   * gmtgrdio_get_image_data (struct GMT_IMAGE
 
 /*! gmt_M_grd_get_size computes grid size including the padding, and doubles it if complex values */
 GMT_LOCAL size_t gmtgrdio_grd_get_size (struct GMT_GRID_HEADER *h) {
-	return ((((h->complex_mode & GMT_GRID_IS_COMPLEX_MASK) > 0) + 1ULL) * h->mx * h->my);
+	return (((((h->complex_mode & GMT_GRID_IS_COMPLEX_MASK) > 0) + 1ULL) * h->mx * h->my) * MAX (h->n_layers, 1));
 }
 
 GMT_LOCAL int gmtgrdio_grd_layout (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, unsigned int complex_mode, unsigned int direction) {
@@ -1953,6 +1953,7 @@ void gmt_set_grddim (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *h) {
 	/* Assumes pad is set and then computes n_columns, n_rows, mx, my, nm, size, xy_off based on w/e/s/n.  */
 	h->n_columns = gmt_M_grd_get_nx (GMT, h);		/* Set n_columns, n_rows based on w/e/s/n and offset */
 	h->n_rows = gmt_M_grd_get_ny (GMT, h);
+	if (h->n_layers == 0) h->n_layers = 1;		/* Default is one grid layer */
 	h->mx = gmt_M_grd_get_nxpad (h, h->pad);	/* Set mx, my based on h->{n_columns,n_rows} and the current pad */
 	h->my = gmt_M_grd_get_nypad (h, h->pad);
 	h->nm = gmt_M_grd_get_nm (h);		/* Sets the number of actual data items */
@@ -2786,6 +2787,87 @@ int gmt_set_outgrid (struct GMT_CTRL *GMT, char *file, bool separate, unsigned i
 	return (false);
 }
 
+int gmtgrdio_init_grdheader (struct GMT_CTRL *GMT, unsigned int direction, struct GMT_GRID_HEADER *header, struct GMT_OPTION *options,
+                             uint64_t dim[], double wesn[], double inc[], unsigned int registration, unsigned int mode) {
+	/* Convenient way of setting a header struct wesn, inc, and registration, then compute dimensions, etc. */
+	double wesn_dup[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, inc_dup[3] = {0.0, 0.0, 0.0};
+	unsigned int n_layers = 1;
+	char *regtype[2] = {"gridline", "pixel"};
+	struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
+	gmt_M_unused(mode);
+
+	if (registration & GMT_GRID_DEFAULT_REG) registration |= GMT->common.R.registration;	/* Set the default registration */
+	registration = (registration & 1);	/* Knock off any GMT_GRID_DEFAULT_REG bit */
+	if (dim && wesn == NULL && inc == NULL) {	/* Gave dimensions instead, set range and inc (1/1) while considering registration */
+		gmt_M_memset (wesn_dup, GMT->common.R.dimension * 2, double);
+		wesn_dup[XHI] = (double)(dim[GMT_X]);
+		wesn_dup[YHI] = (double)(dim[GMT_Y]);
+		inc_dup[GMT_X] = inc_dup[GMT_Y] = 1.0;
+		if (registration == GMT_GRID_NODE_REG) wesn_dup[XHI] -= 1.0, wesn_dup[YHI] -= 1.0;
+		if (dim[GMT_Z] > 1) n_layers = (unsigned int)dim[GMT_Z];
+		GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Grid/Image dimensions imply w/e/s/n = 0/%g/0/%g, inc = 1/1, %s registration, n_layers = %u\n",
+			wesn_dup[XHI], wesn_dup[YHI], regtype[registration], n_layers);
+	}
+	else {	/* Must infer dimension etc from wesn, inc, registration */
+		if (wesn == NULL) {	/* Must select -R setting */
+			if (!GMT->common.R.active[RSET]) {
+				GMT_Report (GMT->parent, GMT_MSG_ERROR, "No w/e/s/n given and no -R in effect.  Cannot initialize new grid\n");
+				GMT_exit (GMT, GMT_ARG_IS_NULL); return GMT_ARG_IS_NULL;
+			}
+		}
+		else	/* In case user is passing header->wesn etc we must save them first as gmt_grd_init will clobber them */
+			gmt_M_memcpy (wesn_dup, wesn, GMT->common.R.dimension * 2, double);
+		if (inc == NULL) {	/* Must select -I setting */
+			if (!GMT->common.R.active[ISET]) {
+				GMT_Report (GMT->parent, GMT_MSG_ERROR, "No increment given and no -I in effect.  Cannot initialize new grid\n");
+				GMT_exit (GMT, GMT_ARG_IS_NULL); return GMT_ARG_IS_NULL;
+			}
+		}
+		else	/* In case user is passing header->inc etc we must save them first as gmt_grd_init will clobber them */
+			gmt_M_memcpy (inc_dup, inc, GMT->common.R.dimension, double);
+		if (dim && dim[GMT_Z] > 1) n_layers = (unsigned int)dim[GMT_Z];
+		if (n_layers == 1 && inc && inc[GMT_Z] > 0.0) {	/* Compute number of layers from equidistant z-range */
+			n_layers = gmt_M_get_n (API->GMT, wesn[ZLO], wesn[ZHI], inc[GMT_Z], 0.0);
+		}
+		if (inc != NULL) {
+			GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Grid/Image dimensions imply w/e/s/n = %g/%g/%g/%g, inc = %g/%g, %s registration, n_layers = %u\n",
+			            wesn_dup[XLO], wesn_dup[XHI], wesn_dup[YLO], wesn_dup[YHI], inc[GMT_X], inc[GMT_Y], regtype[registration], n_layers);
+		}
+	}
+	/* Clobber header and reset */
+	gmt_grd_init (GMT, header, options, false);	/* This is for new grids only so update is always false */
+	if (dim == NULL && wesn == NULL)
+		gmt_M_memcpy (header->wesn, GMT->common.R.wesn, 4, double);
+	else
+		gmt_M_memcpy (header->wesn, wesn_dup, 4, double);
+	if (dim == NULL && inc == NULL)
+		gmt_M_memcpy (header->inc, GMT->common.R.inc, 2, double);
+	else
+		gmt_M_memcpy (header->inc, inc_dup, 2, double);
+	header->registration = registration;
+	/* Copy row-order from R.row_order, if set */
+	if (GMT->common.R.row_order) HH->row_order = GMT->common.R.row_order;
+	/* Mode may contain complex mode information */
+	header->complex_mode = (mode & GMT_GRID_IS_COMPLEX_MASK);
+	HH->grdtype = gmtlib_get_grdtype (GMT, direction, header);
+	gmt_RI_prepare (GMT, header);	/* Ensure -R -I consistency and set n_columns, n_rows in case of meter units etc. */
+	gmt_M_err_pass (GMT, gmt_grd_RI_verify (GMT, header, 1), "");
+	gmt_M_grd_setpad (GMT, header, GMT->current.io.pad);	/* Assign default GMT pad */
+	if (dim) header->n_bands = n_layers;
+	if (GMT->common.R.dimension == 3 && GMT->common.R.inc[GMT_Z] > 0.0) {	/* Compute number of layers from equidistant z-range */
+		if (wesn == NULL)
+			n_layers = gmt_M_get_n (API->GMT, GMT->common.R.wesn[ZLO], GMT->common.R.wesn[ZHI], GMT->common.R.inc[GMT_Z], 0.0);
+		else
+			n_layers = gmt_M_get_n (API->GMT, wesn[ZLO], wesn[ZHI], GMT->common.R.inc[GMT_Z], 0.0);
+		header->n_layers = n_layers;
+	}
+	gmt_set_grddim (GMT, header);	/* Set all dimensions before returning */
+	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) */
+	return (GMT_NOERROR);
+}
+
 int gmt_change_grdreg (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, unsigned int registration) {
 	unsigned int old_registration;
 	double F;

src/gmt_init.c

@@ -8144,6 +8144,7 @@ int gmt_parse_R_option (struct GMT_CTRL *GMT, char *arg) {
 
 	if (GMT->current.setting.run_mode == GMT_MODERN) {	/* Must handle any internal history regarding increments and registration */
 		/* Here, item may be of the form <region>[+I<incs>][+G[P|G][B|T]] if -R was given to a non-plotting module */
+		GMT->common.R.dimension = 2;	/* Since we know it is a grid */
 		if ((c = strstr (item, "+G")) != NULL) {	/* Got grid registration */
 			GMT->common.R.active[GSET] = true;
 			GMT->common.R.registration = strchr(&c[2], 'P') != NULL || strchr(&c[2], 'p') != NULL ? GMT_GRID_PIXEL_REG : GMT_GRID_NODE_REG;
@@ -8225,6 +8226,7 @@ int gmt_parse_R_option (struct GMT_CTRL *GMT, char *arg) {
 		GMT->common.R.wesn[YLO] = orig[GMT_Y] - 0.5 * part * ydim;
 		GMT->common.R.wesn[XHI] = GMT->common.R.wesn[XLO] + xdim;
 		GMT->common.R.wesn[YHI] = GMT->common.R.wesn[YLO] + ydim;
+		GMT->common.R.dimension = 2;	/* Since we know it is a grid */
 		return (GMT_NOERROR);
 	}
 	if ((item[0] == 'g' || item[0] == 'd') && item[1] == '\0') {	/* Check -Rd|g separately in case user has files called d or g */
@@ -8238,6 +8240,7 @@ int gmt_parse_R_option (struct GMT_CTRL *GMT, char *arg) {
 		}
 		GMT->common.R.wesn[YLO] = -90.0;	GMT->common.R.wesn[YHI] = +90.0;
 		gmt_set_geographic (GMT, GMT_IN);
+		GMT->common.R.dimension = 2;	/* Since we know it is a grid */
 		return (GMT_NOERROR);
 	}
 	ptr = item;	/* To avoid compiler warning that item cannot be NULL */
@@ -8324,6 +8327,7 @@ int gmt_parse_R_option (struct GMT_CTRL *GMT, char *arg) {
 		else
 			GMT->current.io.geo.range = GMT_IS_0_TO_P360_RANGE;
 		gmt_set_geographic (GMT, GMT_IN);
+		GMT->common.R.dimension = 2;	/* Since we know it is a grid */
 		return (GMT_NOERROR);
 	}
 	else if (strchr (GMT_LEN_UNITS2, item[0])) {	/* Obsolete: Specified min/max in projected distance units */
@@ -8473,6 +8477,7 @@ int gmt_parse_R_option (struct GMT_CTRL *GMT, char *arg) {
 		GMT_Report (GMT->parent, GMT_MSG_WARNING, "-R with six parameters but no -Jz|Z given - ignore zmin/zmax\n");
 		GMT->common.R.wesn[ZLO] = GMT->common.R.wesn[ZHI] = 0.0;
 	}
+	GMT->common.R.dimension = i / 2;	/* So 2 or 3 */
 	return (error);
 }
 
@@ -14349,6 +14354,7 @@ struct GMT_CTRL *gmt_init_module (struct GMTAPI_CTRL *API, const char *lib_name,
 		GMT_Report (API, GMT_MSG_DEBUG, "Revised options: %s\n", string);
 		GMT_Destroy_Cmd (API, &string);
 	}
+
 	/* Here we can call the rest of the initialization */
 
 	return (gmtinit_begin_module_sub (API, lib_name, mod_name, Ccopy));
@@ -16857,6 +16863,8 @@ struct GMT_CTRL *gmt_begin (struct GMTAPI_CTRL *API, const char *session, unsign
 
 	gmtinit_set_today (GMT);	/* Determine today's rata die value */
 
+	GMT->common.R.dimension = 2;	/* Most likely scenario but will be reset by -R for those cases when it is not a grid */
+
 	return (GMT);
 }
 

src/gmt_resources.h

@@ -385,7 +385,7 @@ enum GMT_enum_grdlen {
 struct GMT_GRID_HEADER {
 	/* Variables we document for the API:
 	 * == Do not change the type of the following three items.
-	 * == They are copied verbatim to the native grid header and must be 4-byte unsigned ints. */
+	 * == They are copied verbatim to any native grid header and must be 4-byte unsigned integers. */
 	uint32_t n_columns;              /* Number of columns */
 	uint32_t n_rows;                 /* Number of rows */
 	uint32_t registration;           /* GMT_GRID_NODE_REG (0) for node grids, GMT_GRID_PIXEL_REG (1) for pixel grids */
@@ -399,7 +399,7 @@ struct GMT_GRID_HEADER {
 	double z_min;                          /* Minimum z value */
 	double z_max;                          /* Maximum z value */
 	double inc[2];                         /* x and y increment */
-	double z_scale_factor;                 /* grd values must be multiplied by this */
+	double z_scale_factor;                 /* grid values must be multiplied by this */
 	double z_add_offset;                   /* After scaling, add this */
 	char x_units[GMT_GRID_UNIT_LEN80];     /* units in x-direction */
 	char y_units[GMT_GRID_UNIT_LEN80];     /* units in y-direction */
@@ -409,21 +409,25 @@ struct GMT_GRID_HEADER {
 	char remark[GMT_GRID_REMARK_LEN160];   /* comments re this data set */
 
 	/* Items not stored in the data file for grids but explicitly used in macros computing node numbers */
-	size_t nm;                  /* Number of data items in this grid (n_columns * n_rows) [padding is excluded] */
-	size_t size;                /* Actual number of items (not bytes) required to hold this grid (= mx * my), per band (for images) */
-	unsigned int bits;          /* Bits per data value (e.g., 32 for ints/floats; 8 for bytes) */
-	unsigned int complex_mode;  /* 0 = normal, GMT_GRID_IS_COMPLEX_REAL = real part of complex grid, GMT_GRID_IS_COMPLEX_IMAG = imag part of complex grid */
-	unsigned int type;          /* Grid format */
-	unsigned int n_bands;       /* Number of bands [1]. Used with IMAGE containers and macros to get ij index from row,col, band */
-	unsigned int mx, my;        /* Actual dimensions of the grid in memory, allowing for the padding */
-	unsigned int pad[4];        /* Padding on west, east, south, north sides [2,2,2,2] */
-	char   mem_layout[4];       /* Three or Four char codes T|B R|C S|R|S (grd) or B|L|P + A|a (img) describing array layout in mem and interleaving */
-	gmt_grdfloat  nan_value;    /* Missing value as stored in grid file */
-	double xy_off;              /* 0.0 (registration == GMT_GRID_NODE_REG) or 0.5 ( == GMT_GRID_PIXEL_REG) */
-	char *ProjRefPROJ4;         /* To store a referencing system string in PROJ.4 format */
-	char *ProjRefWKT;           /* To store a referencing system string in WKT format */
-	int ProjRefEPSG;            /* To store a referencing system EPSG code */
-	void *hidden;               /* Lower-level information for GMT use only */
+	uint32_t n_layers;                     /* Number of layers in a 3-D data cube [1] */
+	size_t nm;                             /* Number of data items in this grid (n_columns * n_rows) [padding is excluded] */
+	size_t size;                           /* Actual number of items (not bytes) required to hold this grid (= mx * my) */
+	unsigned int bits;                     /* Bits per data value (e.g., 32 for ints/floats; 8 for bytes) */
+	unsigned int complex_mode;             /* 0 = normal, GMT_GRID_IS_COMPLEX_REAL = real part of complex grid, GMT_GRID_IS_COMPLEX_IMAG = imag part of complex grid */
+	unsigned int type;                     /* Grid format */
+	unsigned int l_type;                   /* Layer dimension data type */
+	unsigned int n_bands;                  /* Number of bands [1]. Used with IMAGE containers and macros to get ij index from row,col, band */
+	unsigned int mx, my;                   /* Actual dimensions of the grid in memory, allowing for the padding */
+	unsigned int pad[4];                   /* Padding on west, east, south, north sides [2,2,2,2] */
+	char l_units[GMT_GRID_UNIT_LEN80];     /* Name/units for the layer dimension */
+	char mem_layout[4];                    /* Three or Four char codes T|B R|C S|R|S (grd) or B|L|P + A|a (img) describing array layout in mem and interleaving */
+	gmt_grdfloat  nan_value;               /* Missing value as stored in grid file */
+	double xy_off;                         /* 0.0 (registration == GMT_GRID_NODE_REG) or 0.5 ( == GMT_GRID_PIXEL_REG) */
+	double *layer;                         /* Array with layer levels for 3-D cubes [NULL] */
+	char *ProjRefPROJ4;                    /* To store a referencing system string in PROJ.4 format */
+	char *ProjRefWKT;                      /* To store a referencing system string in WKT format */
+	int ProjRefEPSG;                       /* To store a referencing system EPSG code */
+	void *hidden;                          /* Lower-level information for GMT use only */
 };
 
 /* grd is stored in rows going from west (xmin) to east (xmax)

src/gmt_support.c

@@ -6870,7 +6870,7 @@ bool gmt_getinc (struct GMT_CTRL *GMT, char *line, double inc[]) {
 
 	if (!line) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "No argument given to gmt_getinc\n"); return (true); }
 
-	if ((n = gmt_getincn (GMT, line, inc, 2)) < 0) return true;
+	n = gmt_getincn (GMT, line, inc, GMT->common.R.dimension);
 	if (n == 1) {	/* Must copy y info from x */
 		inc[GMT_Y] = inc[GMT_X];
 		GMT->current.io.inc_code[GMT_Y] = GMT->current.io.inc_code[GMT_X];	/* Use exact inc codes for both x and y */