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

[joa-quim]

To start with I don't even understand how the bottom image in #4003 even works. Whatever makes it work must happen in the projection operations because linear projections fail as well. E.g, this fails.

gmt grdimage @earth_day_10m -Baf -R-85/-54/9/26 -JQ10 -P > lixo.ps

But investigating this resurfaced the fact that when reading sub-regions of an image in fact loads the entire image and this is a NO NO.

This PR makes several changes to address that but still doesn't fully work. The PAD is a nightmare. Currently, this works

gmt grdimage @earth_day_10m -Baf -R-85/-54/9/26 -JQ10 -P > lixo.ps

but this not

grdimage @earth_day_10m -Baf -R-85/-54/9/26 -JM10 -P > lixo.ps
grdimage [WARNING]: The image memory layout (TRP ) is of a wrong type. It should be BRPa.
grdimage [ERROR]: gmt_img_project: Input image does not have sufficient (2) padding

If I set the pad = 2, than no complains but image is screwed.

Curiously, this works and no complains with pad = 0.

gmt grdimage @earth_day_10m -Baf -R-85/-54/9/26 -JQ10 -Ilixo_int.grd -P > lixo.ps

[PaulWessel]

FYI, I have a full day of Chair work so maybe afternoon

[PaulWessel]

I wonder if we should first split grdimage into three separate sub-functions:

1. Input data is z-grid only (local or remote), with or without intensity grid or autoderive.
2. input data is red, green blue grids with or without intensity grid [Cannot do autoderive in this case]
3. Input data is an image, with or without intensity grid [cannot do autoderive in this case]

As it is, it is so hard to try to follow the logic in grdimage. While this will set us back some days, might it not be a better and long-overdue update? I have 2 hours now before meetings start up again...

[joa-quim]

I don't know if anyone is still using the 3 r,g,b grids. That's a so old tech.
Breaking in pieces is a so needed step for grdimage ... and others.

[PaulWessel]

With grdmix we have a module that can take 3 grids and write an image. So we could decide to stop accepting 3 grids in grdimage (and grdview) and direct anyone with that type of data to run grdmix first. It would simplify grdimage a bit. I do not think there are tons of people doing 3 grids, and they could then do it that way instead. Alternatively, to not break backwards compatibility, we could have grdimage/view call grdmix first and get that image back.

[joa-quim]

Yep, redirecting to grdmix is a good alternative. It allows us to drop that with no further effort.

[PaulWessel]

So OK for you if I start some hacking later today on grdimage?

[joa-quim]

Sure.

[seisman]

With grdmix we have a module that can take 3 grids and write an image. So we could decide to stop accepting 3 grids in grdimage (and grdview) and direct anyone with that type of data to run grdmix first. It would simplify grdimage a bit. I do not think there are tons of people doing 3 grids, and they could then do it that way instead. Alternatively, to not break backwards compatibility, we could have grdimage/view call grdmix first and get that image back.

If it means a lot of code changes, perhaps we should leave it to 6.2?

[PaulWessel]

Given there will always be bugs, we will never release 6.1.1 unless we just do it. I have not started on grdimage because of work (first week of classes has lots of meetings), so perhaps we should just wrap 6.1.1 as is.

[seisman]

so perhaps we should just wrap 6.1.1 as is.

Yes to me.

[joa-quim]

The large change could be left to after 6.1.1 but we still have #4003 that I've not been able to fix.

And, you'll forgive me but why do we need a 6.1.1 and leave gmtbinstats hanging on the wall for another ~6 months till 6.2? Yes, because if we release 6.1.1 we won't release 6.2 within 2 months. That makes no sense.

[PaulWessel]

I think one reason for 6.1.1 is that it is likely to live in ubuntu and others for a long time before they move to 6.2. So it provides stability. As we are making it easier for people to install (note the auto-download for gshhg and dcw now) we could also add a plug on our website explaining why building the latest has lots of benefits (more bug fixes, more new features) than sticking with a point release.

[PaulWessel]

We also know from past experience that when we release a new module that only I have tried 4-5 times it is likely to have more bugs that your typical module, so having people play with it in master for some time is not a bad idea.

[joa-quim]

Being 6.1.1 or 6.2 makes no difference to Ubuntu. They update versions based on their (slow) release schedule and use whatever is available then. Also, the building argument does not apply to > 50% users (Windowers + many).
No problems in finding new bugs in new modules.

[PaulWessel]

But investigating this resurfaced the fact that when reading sub-regions of an image in fact loads the entire image and this is a NO NO.

Seems like the GMT_Read_Data for images works differently for images and grids. it comes down to gmtapi_import_image versus gmtapi_import_grid. If I recall, I think both correctly read the header (this will return the w/e/s/n for the entire file) but when a grid is read we pass s_obj->wesn to limit the subset, whereas in gmtapi_import_image we do not. Hence the whole thing is read.

[PaulWessel]

I know you dont like the bug-release plan, but we in fact decided this a month ago to do a 6.1.1 and have spent lots of time prepping for that.

[joa-quim]

That was my work in this PR but I'm stuck with the PAD (and more shits with pix-grid-reg probably).

[PaulWessel]

If you make a branch I might be able to look and help as well. More meetings coming up today though.

[seisman]

I think we decided to release major versions every 6 months, and release minor versions if there are any critical issues.

Following that release plan, we're still going to release 6.2.0 after six months of the 6.1.0 release. I don't remember the exact release plan, is it November? During any two major versions, we still release minor versions (like 6.1.1).

[joa-quim]

If you make a branch I might be able to look and help as well. More meetings coming up today though.

This PR is in a branch

[joa-quim]

If we have nothing new, we are not going to release a 6.x.0
If we have new things it's a waste of efforts to leaving them out waiting for a rule that we impose to ourselves.
Point releases with no big bugs forcing it are a lot of work for very little gain.

[PaulWessel]

Well, there has never been a time without new stuff. We keep adding features all the time, and often modules.
I am OK with us revisiting the release policy at the next meeting, but we did agree at that time (at least the majority) so we are sticking with it. To be revisited.

[PaulWessel]

I will try to look at this later today, to many student conferences during the first week of classes...

[PaulWessel]

One issue may be that some of the helper functions, like gmtgrdio_padspace, assumes that the reading function knows how to handle periodic data. Say your image is 0/360 and you want -30/45. Well, for grids we pass in -30/45 and the gmt_nc is smart enough to read the end of a row then padding with the beginning of the row to bet -30/45. I am not sure if gmt_gdalread works the same way? Also, because of padding, the 0/360 request can become -2/362 and it works fine with grids. So perhaps that is why you read the entire image and let the projection deal with it on the GMT side.

[joa-quim]

You're right for the -30/45 case but my example case uses -85/-54 and in gmtlib_read_image() the gmtgrdio_padspace() reads the pad zone from data, but when projecting it still complains that pad is != 2

[PaulWessel]

You can experiment with commenting out the pad check at top or gmt_img_project in gmt_map.c Sorry, have doctors apt and will be gone for part of the day. Paul Wessel, Professor and Chair Dept. of Earth Sciences SOEST, U of Hawaii at Manoa On August 26, 2020 at 4:13:01 AM, Joaquim ([email protected]) wrote: You're right for the -30/45 case but my example case uses -85/-54 and in gmtlib_read_image() the gmtgrdio_padspace() reads the pad zone from data, but when projecting it still complains that pad is != 2 — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#4044 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AGJ7IXZYVFQCLDKH7XS55HTSCUJ6ZANCNFSM4QJ3FDGA> .

[seisman]

We fixed many API bugs after the 6.1.0 release. These fixes are very important for PyGMT. There are many pending PyGMT PRs/features waiting for a new release. Perhaps we can leave the issue #4003 to 6.2 and release GMT 6.1.1 ASAP if we can't fix it in a few days?

[PaulWessel]

Yes, I am OK with that. Hard to predict the timing on #4003. Let's see what @joa-quim says tomorrow and we will go with that.

[joa-quim]

Yes, agree that this issue should not hold the release.

[PaulWessel]

I wish to fix this problem now in 6.2. Currently, grdimage is too complicated to easily debug. Some of the complications are:

• Input: Either one grid, three grids (r,g,b), or one image. These may be remote files or tiles, memory files, or actual files
• Intensity: Optional, but is either a grid, memory reference, constant, or a request to derive from input
• Project: Unless rectangular graticules in the projection, we must warp these items using the projection
• Render: Either we plot PostScript or we do an image directly, writing via GDAL
• Bits: We may convert the input to a 24-bit, 8-bit, or 1-bit image [with -A this becomes 8bit with 2 colors in GDAL]

I think it is best to approach this mess is steps. The one vs three grids is a legacy from before we could use GDAL to read images. Having to deal with three vs one grid complicates our workflow. Because we need to be backwards compatible, I suggest we do the following first:

Detect if we have one or three grids. If three, we call grdmix to obtain the equivalent image, and proceed to plot this image instead.

I can make a test for this, implement the changes, and make sure everything else works after ripping out the 1-vs-3 part both in the code and documentation (I don't think we have any tests for this but we do for grdmix).

Seems and OK start, @joa-quim ?

[PaulWessel]

Hm, I think this PR is being replaced by #4215, no?

[joa-quim]

Most of the changes in this PR are outside grdimage. The only one there may raise a conflict but the rest is probably to maintain.

[PaulWessel]

OK, shall we remove WIP and merge this of is there more testing to be done. You are right, it is barely toughing grdimage.c so should be easy to merge.

[PaulWessel]

After merging in master into this PR I get many failures that all seem to happen in gdal:

The following tests FAILED:
93 - doc/scripts/GMT_images.sh (Failed)
185 - doc/examples/ex16/ex16.sh (Failed)
188 - doc/examples/ex19/ex19.sh (Failed)
220 - doc/examples/ex52/ex52.sh (Failed)
439 - test/grdimage/marbles.sh (Failed)
442 - test/grdimage/png_image.sh (Failed)
443 - test/grdimage/readwrite_withgdal.sh (Failed)
448 - test/grdimage/url_map.sh (Failed)
667 - test/psimage/psimage.sh (Failed)
739 - test/pstext/filled_font.sh (Failed)
845 - test/psxy/tiling.sh (Failed)

For instance, ex19 crashes on middle panel when we use coast to fill using a pattern image. Crash happens in gmt_gdalread.c line 1058 because nYSize = nXSize == 0.

	/* 16 Mb worth of rows */
	nRowsPerBlock = MIN(nYSize, (int)(1024 * 1024 * 16 / (nXSize * nPixelSize)));
	nBlocks = (int)ceil((float)nYSize / nRowsPerBlock);

You can recreate the crash using this branch and running

gmt coast -Dc [email protected]+r100 -Rd -JI0/16c -pdf map

so you can see what is happening.

[PaulWessel]

Hi @joa-quim, perhaps merge in master into this branch and give it a check again. It has lingered as a WIP for a while now.

[joa-quim]

I have some urgent bugs to fix in Mirone. After it I'll give this another try.

[joa-quim]

Still

grdimage @earth_day_10m -Baf -R-85/-54/9/26 -JM10 -P > lixo.ps
grdimage [ERROR]: gmt_img_project: Input image does not have sufficient (2) padding

[joa-quim]

Don't know what to do to fix this. Spend many hours debugging but just can't get rid of the pad shit. If the -R fits all inside then it's working already. One would only have to find a way to tell gmt_img_project() that the pad is there in fact. But when the -R touches one or more of the data edges then a true pad is need but gdalread doesn't allow reading with a -R that exceeds data domain (that -R is set by gmtgrdio_padspace()).

[PaulWessel]

I think for grids without pad, the equivalent gmt_grd_project makes a temp grid with the pad so the interpolation can happen. Perhaps the same can be done for the images?

[joa-quim]

Duplicating the entire image because of the pad is too crazy wasteful.

[PaulWessel]

Once I have time to come up for air I will debug to re-educate myself on what is happening there.

[Esteban82]

It would be possible that this was fixed with #5798 ?

[PaulWessel]

Yes, I think so. Basically the branch I worked on was meant to replace this branch which probably will be deleted. Let me know if that I not correct, @joa-quim. I looked at it and copied a few things into the other branch now merged into master as well as #5819.

[PaulWessel]

Hi @joa-quim I am pretty sure this branch and PR can be deleted, no? This was dealt with a month or more ago (image subsets).

[joa-quim]

Yes,to close.