Implement image subsets within GMT and via grdcut #5798

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged

PaulWessel merged 34 commits into master from grdcut-image

Sep 28, 2021

doc/rst/source/grdcut.rst

-Original file line number
+Diff line change
@@ Expand Up @@
     indirectly via a range check on the node values or via distances from a
     given point. Finally, you can use **-J** for oblique projections to determine
     the corresponding rectangular **-R** setting that will give a grid that fully
-    covers the oblique domain.
+    covers the oblique domain.  **Note**: If input grid is actually an image, then
+    only options **-R** and **-G** are supported, i.e., you can cut out a sub-region only.
     Complementary to **grdcut** there is :doc:`grdpaste`, which
     will join together two grid files along a common edge.
@@ Expand Down @@

src/gmt_api.c

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -4368,17 +4368,15 @@ GMT_LOCAL bool gmtapi_expand_index_image (struct GMT_CTRL *GMT, struct GMT_IMAGE
  
    		for (c = 0; c < 3; c++) off[c] = c * h->size;

    		for (node = 0; node < h->size; node++) {	/* For all pixels, including the pad */

    			index = I->data[node];	/* Pixel index into color table */

    			start_c = index * 4;	/* Start of the r,g,b entry in the color table for this index */

    			for (c = 0; c < 3; c++) data[node+off[c]] = I->colormap[start_c+c];	/* Place r,g,b in separate bands */

    			for (c = 0; c < 3; c++) data[node+off[c]] = gmt_M_get_rgba (I->colormap, index, c, I->n_indexed_colors);	/* Place r,g,b in separate bands */

    		}

    	}

    	else {	/* Pixel interleave */

    		uint64_t k;

    		strncpy (h->mem_layout, "TRP ", 4);	/* Fill out red, green, and blue pixels */

    		for (node = k = 0; node < h->size; node++) {	/* For all pixels, including the pad */

    			index = I->data[node];	/* Pixel index into color table */

    			start_c = index * 4;	/* Start of the r,g,b entry in the color table for this index */

    			for (c = 0; c < 3; c++, k++) data[k] = I->colormap[start_c+c];	/* Place r,g,b in separate bands */

    			for (c = 0; c < 3; c++, k++) data[k] = gmt_M_get_rgba (I->colormap, index, c, I->n_indexed_colors);	/* Place r,g,b in separate bands */

    		}

    		/* If neither TRB or TRP we call for a changed layout, which may or may not have been implemented */

    		GMT_Change_Layout (GMT->parent, GMT_IS_IMAGE, GMT->parent->GMT->current.gdal_read_in.O.mem_layout, 0, I, NULL, NULL);

    @@ -4400,7 +4398,7 @@ int gmtlib_ind2rgb (struct GMT_CTRL *GMT, struct GMT_IMAGE **I_in) {
  
    	   called by gmtapi_import_image, there are other cases when we need also to convert from indexed to RGB.

    	   For example in grdimage when the image was sent in via an external wrapper. In this case the code flow goes

    	   through gmtapi_get_image_data() (in GMT_Read_Data -> gmtapi_pass_object (API, S_obj, family, mode, wesn))

    	   and deliver that Image object directly to the calling module amd may thus have indexed pixels.

    	   and deliver that Image object directly to the calling module and may thus have indexed pixels.

    	*/

    	struct GMT_IMAGE* Irgb = NULL;

    	if ((*I_in)->header->n_bands == 1 && (*I_in)->n_indexed_colors > 0) {		/* Indexed image, convert to RGB */

    @@ -4428,7 +4426,7 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob
  
    	 */

    	int item, new_item, new_ID;

    	bool done = true, via = false, must_be_image = true, no_index = false;

    	bool done = true, via = false, must_be_image = true, no_index = false, bc_not_set = true;

    	uint64_t i0, i1, j0, j1, ij, ij_orig, row, col;

    	unsigned int both_set = (GMT_CONTAINER_ONLY | GMT_DATA_ONLY);

    	double dx, dy, d;

    @@ -4490,20 +4488,24 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob
  
    			/* 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 (GMT->common.R.active[RSET] && !S_obj->region) { /* subregion not passed to object yet */

                    gmt_M_memcpy (S_obj->wesn, GMT->common.R.wesn, 4U, double);

                    S_obj->region = true;

                }

    			size = gmtapi_set_grdarray_size (GMT, I_obj->header, mode, S_obj->wesn);    /* Get array dimension only, which includes padding. DANGER DANGER JL*/

    			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->type <= GMT_UCHAR)

    					I_obj->data = gmt_M_memory (GMT, NULL, I_obj->header->size * I_obj->header->n_bands, unsigned char);

    					I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned char);

    				else if (I_obj->type <= GMT_USHORT)

    					I_obj->data = gmt_M_memory (GMT, NULL, I_obj->header->size * I_obj->header->n_bands, unsigned short);

    					I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned short);

    				else if (I_obj->type <= GMT_UINT)

    					I_obj->data = gmt_M_memory (GMT, NULL, I_obj->header->size * I_obj->header->n_bands, unsigned int);

    					I_obj->data = gmt_M_memory (GMT, NULL, size * I_obj->header->n_bands, unsigned int);

    				else {

    					GMT_Report (API, GMT_MSG_ERROR, "Unsupported image data type %d\n", I_obj->type);

    					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*/

    				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);

    @@ -4512,8 +4514,11 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob
  
    			else if (gmt_M_err_pass (GMT, gmtlib_read_image (GMT, S_obj->filename, I_obj, S_obj->wesn,

    				I_obj->header->pad, mode), S_obj->filename))

    				return_null (API, GMT_IMAGE_READ_ERROR);

    			if (gmt_M_err_pass (GMT, gmtlib_image_BC_set (GMT, I_obj), S_obj->filename))

    				return_null (API, GMT_IMAGE_BC_ERROR);	/* Set boundary conditions */

    			if (I_obj->n_indexed_colors == 0) {	/* May set the BCs */

    				if (gmt_M_err_pass (GMT, gmtlib_image_BC_set (GMT, I_obj), S_obj->filename))

    					return_null (API, GMT_IMAGE_BC_ERROR);	/* Set boundary conditions */

    				bc_not_set = false;

    			}

    			IH = gmt_get_I_hidden (I_obj);

    			IH->alloc_mode = GMT_ALLOC_INTERNALLY;

    #else

    @@ -4695,11 +4700,16 @@ GMT_LOCAL struct GMT_IMAGE *gmtapi_import_image (struct GMTAPI_CTRL *API, int ob
  
    	if (done) S_obj->status = GMT_IS_USED;	/* Mark as read (unless we just got the header) */

    #ifdef HAVE_GDAL

    	if (no_index && gmtapi_expand_index_image (API->GMT, I_obj, &Irgb)) {	/* true if we have a read-only indexed image and we had to allocate a new one */

    		if (GMT_Destroy_Data (API, &I_obj) != GMT_NOERROR) {

    			return_null (API, API->error);

    	if (no_index) {   /* true if we have an indexed image and we had to allocate a new one */

    		if (gmtapi_expand_index_image (API->GMT, I_obj, &Irgb)) {   /* true if we have a read-only indexed image and we had to allocate a new one */

    			if (GMT_Destroy_Data (API, &I_obj) != GMT_NOERROR) {

    				return_null (API, API->error);

    			}

    			I_obj = Irgb;

    		}

    		I_obj = Irgb;

    		/* If we were unable to set BCs earlier we must do it now */

    		if (bc_not_set && gmt_M_err_pass (GMT, gmtlib_image_BC_set (GMT, I_obj), S_obj->filename))

    			return_null (API, GMT_IMAGE_BC_ERROR);	/* Failed to set boundary conditions */

    	}

    #endif

    @@ -12849,10 +12859,12 @@ 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 */

        else if (!strncmp (module, "grdcut", 6U) && GMT_Find_Option (API, 'D', *head)) {

             deactivate_output = true;   /* Turn off implicit output since none is in effect, only secondary -D output */

        }

        /* 1x. Check if grdcut is just getting information and if input is a grid or image */

        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 */

       }

    	/* 2a. Get the option key array for this module */

    	key = gmtapi_process_keys (API, keys, type, *head, n_per_family, &n_keys);	/* This is the array of keys for this module, e.g., "<D{,GG},..." */

src/gmt_customio.c

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -1909,15 +1909,15 @@ int gmt_gdal_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt
  
    			/* Here we assume that all pad[0] ... pad[3] are equal. Otherwise ... */

    			to_gdalread->mini_hdr.active = false;	/* Undo above setting */

    			to_gdalread->p.active = true;

    			to_gdalread->p.pad = (int)pad[XLO];

    			gmt_M_memcpy (to_gdalread->p.pad, pad, 4U, unsigned int);

    		}

    		/* OK, now test if we are under the condition of #403 (very small grids).

    		   If yes, undo the mini_hdr solution ... and hope no more troubles arise. */

    		if (to_gdalread->mini_hdr.active && (header->n_columns <= 4 || header->n_rows <= 4)) {

    			to_gdalread->mini_hdr.active = false;

    			to_gdalread->p.active = true;

    			to_gdalread->p.pad = (int)pad[XLO];

    			gmt_M_memcpy (to_gdalread->p.pad, pad, 4U, unsigned int);

    		}

    	}

    	if (HH->pocket) {	/* Have a band request. */

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement image subsets within GMT and via grdcut #5798

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!