WIP Try to adress #4003 as well as more underlying issues

src/gmt_api.c

@@ -4340,7 +4340,7 @@ GMT_LOCAL struct GMT_IMAGE * gmtapi_import_image (struct GMTAPI_CTRL *API, int o
 			/* Here we will read the image data themselves. */
 			/* To get a subset we use wesn that is not NULL or contain 0/0/0/0.
 			 * Otherwise we extract the entire file domain */
-			if (!I_obj->data) {	/* Array is not allocated yet, do so now. We only expect header (and possibly w/e/s/n subset) to have been set correctly */
+			if (!I_obj->data && !gmt_M_is_subset (API->GMT, I_obj->header, S_obj->wesn)) {	/* Array is not allocated yet. Allocate now if a subset was not required. If yes allocation is done in gdlaread. We only expect header to have been set correctly */
 				if (I_obj->type <= GMT_UCHAR)
 					I_obj->data = gmt_M_memory (GMT, NULL, I_obj->header->size * I_obj->header->n_bands, unsigned char);
 				else if (I_obj->type <= GMT_USHORT)
@@ -4352,8 +4352,8 @@ GMT_LOCAL struct GMT_IMAGE * gmtapi_import_image (struct GMTAPI_CTRL *API, int o
 					return_null (API, GMT_NOT_A_VALID_TYPE);
 				}
 			}
-			else {	/* Already have allocated space; check that it is enough */
-				size = gmtapi_set_grdarray_size (GMT, I_obj->header, mode, S_obj->wesn);	/* Get array dimension only, which includes padding. DANGER DANGER JL*/
+			else if (I_obj->data) {	/* Already have allocated space; check that it is enough */
+				size = gmtapi_set_grdarray_size (GMT, I_obj->header, mode, S_obj->wesn);	/* Get array dimension only, which includes padding. */
 				if (size > I_obj->header->size) return_null (API, GMT_IMAGE_READ_ERROR);
 			}
 			GMT_Report (API, GMT_MSG_INFORMATION, "Reading image from file %s\n", S_obj->filename);

src/gmt_grdio.c

@@ -3568,24 +3568,30 @@ int gmtlib_read_image (struct GMT_CTRL *GMT, char *file, struct GMT_IMAGE *I, do
 	from_gdalread = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALREAD_OUT_CTRL);
 
 	if (GMT->common.R.active[RSET]) {
-		snprintf (strR, GMT_LEN128, "%.10f/%.10f/%.10f/%.10f", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI],
-		          GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
+		//snprintf (strR, GMT_LEN128, "%.10f/%.10f/%.10f/%.10f", GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI],
+		          //GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
+		snprintf (strR, GMT_LEN128, "%.10f/%.10f/%.10f/%.10f", P.wesn[XLO], P.wesn[XHI], P.wesn[YLO], P.wesn[YHI]);
 		to_gdalread->R.region = strR;
-		/*to_gdalread->R.active = true;*/	/* Wait until we really know how to use it */
+		to_gdalread->registration.val = I->header->registration;	/* Due to pix-reg only by GDAL we need to inform it about our reg type */
+		to_gdalread->registration.x_inc = I->header->inc[GMT_X];
+		to_gdalread->registration.y_inc = I->header->inc[GMT_Y];
+		to_gdalread->R.active = true;	/* Wait until we really know how to use it */
 	}
 
 	if (HH->pocket) {				/* See if we have a band request */
 		to_gdalread->B.active = true;
 		to_gdalread->B.bands = HH->pocket;	/* Band parsing and error testing is done in gmt_gdalread */
 	}
 
-	to_gdalread->p.pad = (int)pad[0];	/* Only 'square' padding allowed */
-	to_gdalread->p.active = (pad[0] > 0);
+	to_gdalread->p.pad = (int)P.pad[0];	/* Only 'square' padding allowed */
+	to_gdalread->p.active = (P.pad[0] > 0);
 	to_gdalread->I.active = true;		/* Means that image in I->data will be BIP interleaved */
 
 	/* Tell gmt_gdalread that we already have the memory allocated and send in the *data pointer */
-	to_gdalread->c_ptr.active = true;
-	to_gdalread->c_ptr.grd = I->data;
+	if (I->data) {		/* Otherwise (subregion) memory is allocated in gdalread */
+		to_gdalread->c_ptr.active = true;
+		to_gdalread->c_ptr.grd = I->data;
+	}
 
 	if (gmt_gdalread (GMT, file, to_gdalread, from_gdalread)) {
 		GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR reading image with gdalread.\n");
@@ -3597,6 +3603,27 @@ int gmtlib_read_image (struct GMT_CTRL *GMT, char *file, struct GMT_IMAGE *I, do
 		return (GMT_GRDIO_READ_FAILED);
 	}
 
+	if (!I->data) {
+		if      (from_gdalread->UInt8.active)  I->data = from_gdalread->UInt8.data;
+		//else if (from_gdalread->UInt16.active) I->data = from_gdalread->UInt16.data;
+		//else if (from_gdalread->Int16.active)  I->data = from_gdalread->Int16.data;
+		//else if (from_gdalread->Int32.active)  I->data = from_gdalread->Int32.data;
+		//else if (from_gdalread->UInt32.active) I->data = from_gdalread->UInt32.data;
+		else {
+			GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR reading image with gdalread.\n");
+			return (GMT_GRDIO_READ_FAILED);
+		}
+		gmt_M_memcpy (I->header->wesn, from_gdalread->hdr, 4, double);
+		I->header->n_columns = from_gdalread->RasterXsize;
+		I->header->n_rows    = from_gdalread->RasterYsize;
+		I->header->nm = gmt_M_get_nm (GMT, I->header->n_columns, I->header->n_rows);
+		// TEMP while the pad shit is not solved
+		pad[0] = pad[1] = pad[2] = pad[3] = 0;
+		I->header->mx   = I->header->n_columns + pad[0] + pad[1];
+		I->header->my   = I->header->n_rows + pad[2] + pad[3];
+		I->header->size = (size_t)I->header->mx * I->header->my;
+	}
+
 	if (to_gdalread->O.mem_layout[0]) 	/* If a different mem_layout request was applied in gmt_gdalread than we must update */
 		gmt_strncpy(I->header->mem_layout, to_gdalread->O.mem_layout, 4);
 

src/grdimage.c

@@ -1253,7 +1253,7 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
 
 		/* Read in the the entire image that is to be mapped */
 		GMT_Report (API, GMT_MSG_INFORMATION, "Allocate memory and read image file %s\n", Ctrl->In.file);
-		if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA | GMT_IMAGE_NO_INDEX, NULL, Ctrl->In.file, NULL)) == NULL) {
+		if ((I = GMT_Read_Data (API, GMT_IS_IMAGE, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA | GMT_IMAGE_NO_INDEX, GMT->common.R.wesn, Ctrl->In.file, NULL)) == NULL) {
 			Return (API->error);
 		}
 		grid_registration = I->header->registration;	/* This is presumably pixel registration since it is an image */