Implement image subsets within GMT and via grdcut

[PaulWessel]

As a way to debug and implement subsets of images, I am letting gmt grdcut -I get the ability to cut a subset out of an image. This also makes it easier to debug this feature which is still not working in GMT (and affects grdimage, for instance, where the entire image must be read in, leading to excessive run times and memory requirements for high-resolution global grids). Derives bits and pieces from stalled PR #4004. Note: I also updated the KEYS and GMT_Encode_Options so it should be possible to receive the image back from the externals as well.

Trying this command:

gmt grdcut @earth_day_01m -R0/10/0/10 -I -Gcut.png

Yields an image with the correct region, increment meta data but likely bad info on the order or RGB, yet the pixels are in the right layout (since Africa outline etc is correct) and the region matches that of a pscoast plot:

Hoping @joa-quim can advise where to make changes in this PR.

Note: This PR allocates memory only for the subset and pad, and currently grdimage commands like the one below will crash in the PR because it assumes a full image. This will be dealt with once we correctly can read subsets of images. However, in master, this comand

gmt grdimage @earth_day_01m -R0/10/0/10 -JQ5.07614213198c -png out

yields this (same-size) image (for comparison with above):

Clearly, this is an incorrect region despite having the RGB in the right order, so this seems to be a bug in master. I suggest we work on getting the subset to work correctly via grdcut and then the rest will be simpler.

[maxrjones]

Not sure that I understand your last post - are you saying to ignore the problems with the colors?

[PaulWessel]

Not sure that I understand your last post - are you saying to ignore the problems with the colors?

No, there are two things:

1. In this PR, looks like grdimage with an image crashes due to the changes we are trying in the i/o for images. This will go away when we are finished.
2. I did not expect master to give the wrong image but I think @federico pointed out some problems with grdimage using the earth_day stuff earlier so maybe the same. I don't really want to stop this PR to fix that in master since the real solution is the one we will get from this PR.

[PaulWessel]

Hi @joa-quim, towards the end of gmtlib_read_image in gmt_grdio.c there is

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

In our case, the image I has TRBa which the to_gdalread version says BRPa, so we copy it over. This gives the above wrong coloring. Should it be another string?

[PaulWessel]

OK, this is an indexed image so not wrong R/G/B order, but upon comparing the image produced with grdcut to the equivalent derived via

gdal_translate -srcwin 10800 4800 600 600 ~/.gmt/server/earth/earth_day/earth_day_01m_p.tif gdal.tif

I find that the two colortables are just very different. How can we get that wrong?

[PaulWessel]

Stepping through to see what happens to the indexed colors.

Line 539 gmt_gdalread.c: We are getting the colors one by one from GDAL and storing in the Ctrl->ColorMap 4x256 array. For each color i = 0, 1, ...you do this

GDALGetColorEntryAsRGB(hTable, i, &sEntry);
Ctrl->ColorMap[j++] = sEntry.c1;
Ctrl->ColorMap[j++] = sEntry.c2;
Ctrl->ColorMap[j++] = sEntry.c3;
Ctrl->ColorMap[j++] = sEntry.c4;

That means the RGBA for the first entry in the colormap is the first 4 values in Ctrl->ColorMap. For my example that is 224,164,84,255.

Then in Line 177 (gmt_gdalwrite.c) you place these values into another array, this time 4x256 float:

for (k = 0; k < n_colors; k++) {
	to_GDALW->C.cpt[k] = I->colormap[k] / 255.0;
	to_GDALW->C.cpt[k + n_colors]   = gmt_M_is255(I->colormap[k + n_colors]);
	to_GDALW->C.cpt[k + n_colors*2] = gmt_M_is255(I->colormap[k + n_colors*2]);
}

So here, let k = 0. Then you get the red (colormap[0]) but not the green from the first (since that would be colormap[1], divide by 255 and store it in the same place in the cpt array. I guess this is OK but not sure why you do it this way since it does not follow the structure of colormap.

Finally, just before we call GDAL to write (line 375) you do this:

for (i = 0; i < nColors; i++) {
	sEntry.c1 = (short)(gmt_M_s255(ptr[i]));
	sEntry.c2 = (short)(gmt_M_s255(ptr[i+nColors]));
	sEntry.c3 = (short)(gmt_M_s255(ptr[i+2*nColors]));
	sEntry.c4 = (short)255;
	GDALSetColorEntry(hColorTable, i, &sEntry);
}

Again, try i = 0. This picks supposedly an r, g, b from the ptr (which is the cpt above) and then loads up the sEntry structure. For i = 0 these become 224,72,72,255. So no longer 224,164,84,255. I think that is why the image looks junk. Please let me know what you were trying to do here, @joa-quim

[joa-quim]

I think there must be missing some step in the middle. While the first case stores the color matrix as row-major the other two are accessing it as column major. The column major ordering is quite likely because when that part was written the main idea was to interface with the MEX.

[PaulWessel]

OK, I can ensure the ColorMap is always stored internally in column order throughout and see how that goes.

[PaulWessel]

There are other places in GMT where the assumtion is that colormap is rows, e.g. in gmtapi_expand_index_image we do

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

I guess I will need to go through all of GMT where ->colormap is referenced and ensure it is set and used as columns. It is not that many places I think.

[PaulWessel]

OK, I have gotten grdcut -I to work. Not tested extensively on different images though, just earth_day, but output matches gdal_translate. This is reading and returning an image with the pad (as expected for instance in grdimage) but when we write we strip off the pad.

gmt grdcut @earth_day_01m -R-5/10/0/10 -I -Gcut.tif

Of course, it cannot handle periodic boundary at ±180:

gmt grdcut @earth_day_01m -R170/190/0/10 -I -Gcut.tif
grdcut [ERROR]: Computed -srcwin falls outside raster size of 21600x10800.
grdcut [ERROR]: ERROR reading image with gdalread
```.

That may require the use of 

`gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_in, &height_in, &first_col, &last_col, &first_row, &last_row, &actual_row), HH->name);`

and reading GDAL twice.  Not sure what to do about this.

[joa-quim]

The read sub-regions part originated long ago in Mirone and either the number of options were more limited, or more likely my knowledge. It seems we should be using -a_ullr instead of -srcwin since former let us specify the coordinates in geogs while latter uses only row/col. GDAL can also deal with round Earth but it often needs to be told to do so.

[maxrjones]

Would it be possible to get this to work without the -I argument being required, similar to how grdimage accepts grid | image without any extra user input?

[maxrjones]

The read sub-regions part originated long ago in Mirone and either the number of options were more limited, or more likely my knowledge. It seems we should be using -a_ullr instead of -srcwin since former let us specify the coordinates in geogs while latter uses only row/col. GDAL can also deal with round Earth but it often needs to be told to do so.

This branch also has a lot of test failures caused by problems similar to 'grdimage [ERROR]: Computed -srcwin falls outside raster size of ...', which could be caused by this.

[PaulWessel]

Would it be possible to get this to work without the -I argument being required, similar to how grdimage accepts grid | image without any extra user input?

Yes, probably. Did -I in this branch just to get things going.

[PaulWessel]

Yes, it looks like we need to upgrade the gdalread to specify geographic coordinates instead when the image is global. Can @joa-quim get us started on that, i.e., how to pass different arguments with coordinates instead?

[Esteban82]

Is this suppose to solve #4003 ? I try the following script. For some resolutions, it worked. But failed for other.
Are you working on this @PaulWessel ?

I try it before the "update grdcut.c" commit.

Failed: 01d, 20m. 06m, 05m, 03m, 02m
Ok: 30m, 15m, 10m, 04m, 01m

RES=01m
gmt begin BlueMarble_Intens png
	gmt grdgradient @earth_relief_$RES -Nt1.00 -A45 -GIntens.nc -R-85/-54/9/26
	gmt grdimage @earth_day_$RES -IIntens.nc -Baf -JM10c 
gmt end show

[PaulWessel]

I see the 01d crashes. Will have a look - I did not test specifically with intensity. We may turn your example into a test once it is working.

[Esteban82]

Now I try it after the commit.

These are the errors (indicated by resolution of the grid/image):

02m, 03m and 01d:
munmap_chunk(): invalid pointer
image.sh: línea 6: 37164 Abortado (`core' generado) gmt grdimage @earth_day_$RES -IIntens.nc -Baf -JM10c

05m and 06m:
free(): invalid next size (normal)
image.sh: línea 6: 37194 Abortado (`core' generado) gmt grdimage @earth_day_$RES -IIntens.nc -Baf -JM10c

20m:
ERROR: Caught signal number 11 (Segmentation fault) at
/lib/x86_64-linux-gnu/libc.so.6(+0x9b941)[0x7f746cc9b941]
[0x0]
Stack backtrace:
/usr/local/lib/libgmt.so.6(sig_handler+0x2ad)[0x7f746ce5475d]
/lib/x86_64-linux-gnu/libc.so.6(+0x46210)[0x7f746cc46210]
/lib/x86_64-linux-gnu/libc.so.6(+0x9b941)[0x7f746cc9b941]
...

[PaulWessel]

This may drag into tomorrow as I have been busy moving my UH home office back to UH today. Setting up my computer rig in the office for the first time since March 2020...

[maxrjones]

The three test differences look fine.

[PaulWessel]

I think I just fixed that @Esteban82 - update and give it a check.

[Esteban82]

I think I just fixed that @Esteban82 - update and give it a check.

Ok. I will check them as soon as posible.

[maxrjones]

I am getting a segfault in a 6.2 script that I have used to plot images created by grdmix. I'll send an email with more details so that I can share the data since it's not using a remote file.

[Esteban82]

I think I just fixed that @Esteban82 - update and give it a check.

It works with all the resolution, even 30s.

Is there any other any test besides Meghan issues.

[PaulWessel]

Thanks, I have added a test based on @Esteban82 example above (same plus one that straddles the Dateline) and it passes. Also fixed the things @meghanrjones mentioned related to her non-global large TIFF. Hopefully we are ready to approve.

[joa-quim]

I'm getting a crash when trying to visualize and indexed image whose that is supposed to have transparency set to one of the rows of the cpt. That is done a bit under the hood by setting

img.n_colors = n_colors * 1000 + Int32(alpha_ind)

However, it also crashes in master so the breaking must have occurred after I implemented it but before the work in this branch. So it should not be a merge stopper. I'll look to what's happening when I have time.