Add grid sampling as new algorithm to grdfill

[PaulWessel]

This PR introduces -Aggrid as a new algorithm, which will sample the given grid at the nodes with NaN in the main input grid. I have added a test that demonstrates it works well. I realized this could be useful in working with David Sandwell on trying to speed up his 5760 x 6912 surface commands in GMTSAR - this algorithm is a partial answer for filling holes. Test output looks like this - also exercises Roman numerals in subplot:

i) original data, ii) mask,iii) masked data, iv) same with image-masking, v) filled grid, vi) differences

[PaulWessel]

OK, we scaled by the test to make it simpler and stored the holed data set in cache, per @meghanrjones suggestion. The server has now synced up the cache so the tests can run again. The plot now looks like this:

DVC has been updated and tests pass. I think this is good to go now.

[seisman]

The grdfill.test fails on Windows:

C:\Program Files\GraphicsMagick-1.3.32-Q8\gm.exe compare: image difference exceeds limit (0.0355743 > 0.003).
test/grdfill/gridfill.ps: RMS Error = 0.0356 [FAIL]
memtrack errors: 0
exit status: 1

Here is the diff image:

Here is the generated PDF file on Windows: gridfill.pdf

[PaulWessel]

The PDF looks reasonable but clearly they differ. Hopeless to know why. It is simply passing x,y arrays to grdtrack and @earth_relief_01d. Funny that one node inside the circle agrees exactly. Maybe @joa-quim could run this test locally and if he gets the same then post the grid that is being made by grdfill so I can compare with what I see.

[seisman]

If it helps you, here is the grid file generated by the CI new.grd.zip

[PaulWessel]

Certainly an interesting one:

gmt grdcut @earth_relief_01d -R0/10/0/10 -Gb.grd
echo 7.5 1.5 | gmt grdtrack -Gb.grd -nn
7.5	1.5	-1169.5
echo 7.5 1.5 | gmt grdtrack -G@earth_relief_01d -R0/10/0/10 -nn
7.5	1.5	-2347.5

Since the script is using the last scheme (passing the remote file to grdfill) there must be something funny here. BTW, the Windows produced grid gives:

echo 7.5 1.5 | gmt grdtrack -Gnew.grd -nn
7.5	1.5	-2701.41601562

while my output form the script gives -1169.5. There is an Island really close to this point so shallow makes sense. I think the grdcut grid is correct so it should be shallow here and for some reason on Windows we are getting another node.

If you or @joa-quim could modify the script to do

gmt grdcut @earth_relief_01d -R0/10/0/10 -Gb.grd
gmt grdfill @earth_relief_20m_holes.grd -Agb.grd -Gnew.grd

then my guess is it will give the same result as me (and pass the test) and then we know that the passing of remote files to grdfill (which passes those args to grdtrack) is suspect.

[PaulWessel]

Hi @joa-quim, might you be able to replace that grdfill command with the two at the end of the previous message and see if that works well for Windows? That way I know what to look for when I debug gmt grdtrack -G@earth_relief_01d -R0/10/0/10 -nn

[joa-quim]

Not sure I got it. Is this what you asked?

gmt grdcut @earth_relief_01d -R0/10/0/10 -Gb.grd
gmt grdfill @earth_relief_20m_holes.grd -Agb.grd -Gnew.grd
echo 7.5 1.5 | gmt grdtrack -Gnew.grd -nn
7.5     1.5     -1169.5

[joa-quim]

And, BTW, the tests pass here.

ctest -R grdfill
Test project C:/v/build
    Start 436: test/grdfill/constfill.sh
1/4 Test #436: test/grdfill/constfill.sh ........   Passed    3.26 sec
    Start 437: test/grdfill/nnfill.sh
2/4 Test #437: test/grdfill/nnfill.sh ...........   Passed    2.33 sec
    Start 438: test/grdfill/showregions.sh
3/4 Test #438: test/grdfill/showregions.sh ......   Passed    2.32 sec
    Start 439: test/grdfill/splinefill.sh
4/4 Test #439: test/grdfill/splinefill.sh .......   Passed    2.64 sec

[PaulWessel]

And what does this give

echo 7.5 1.5 | gmt grdtrack -G@earth_relief_01d -R0/10/0/10 -nn

[seisman]

@PaulWessel For mysterious reasons, the Windows CI now passes for the latest commit (https://github.com/GenericMappingTools/gmt/actions/runs/2300473900).

[PaulWessel]

Funny, so maybe no problem? Of maybe random since I think the things I showed earlier are problematic.

[joa-quim]

Sorry, the test in question was not executed (possibly need to wipe out the build and restart). If I run from within the test dir I get errors but it creates a new.grd grid and with it

echo 7.5 1.5 | gmt grdtrack -Gnew.grd -nn
7.5     1.5     -2701.41601562

So the gridfill.sh test is not generating the same grid as

gmt grdfill @earth_relief_20m_holes.grd -Agb.grd -Gnew.grd

[joa-quim]

... but now I run the test script with GIT bash, no errors reported and the grid is different, correct I assume because

echo 7.5 1.5 | gmt grdtrack -Gnew.grd -nn
7.5     1.5     -1169.5

The error with another shell was

gmt [ERROR]: Not available in classic mode
grdimage [ERROR]: Option -c is not a recognized common option
grdimage [ERROR]: Option -c is not a recognized common option
gmt [ERROR]: Not available in classic mode
gmt [ERROR]: Shared GMT module not found: colorbar

[joa-quim]

echo 7.5 1.5 | gmt grdtrack -G@earth_relief_01d -R0/10/0/10 -nn
7.5     1.5     -2347.5

[PaulWessel]

Thanks, so I need to look at that grdtrack with remote file.

[PaulWessel]

OK, mystery solved. Here is what is going on (not sure of best solution yet):

If registration is not selected then we always pick the pixel registered remote data set when making plots. However, for data processors such as grdtrack we would run into the 1/2 pixel gap near the border when sampling so we switch the default to gridline registration here. Hence

gmt grdtrack -G@earth_relief_01d -R0/10/0/10 -nn t.txt

becomes

gmt grdtrack -G@earth_relief_01d_g -R0/10/0/10 -nn t.txt

whereas

gmt grdcut @earth_relief_01d -R0/10/0/10 -Gb.grd

becomes

gmt grdcut @earth_relief_01d_p -R0/10/0/10 -Gb.grd

Now, none of this is ideal I think. There are good reason for switching to gridline when running grdtrack on a remote grid without specific registration. A question might be if we should do this for all non-plotting grd modules, i.e., grdcut, grd??? etc. At least that would be consistent. Then, this grdfill job would at least select the gridline registered version as will the grdtrack call inside it.

My proposal would be to default to gridline in grd processors and stay with pixel in plotters. The reason is that for plot it does not really matter too much but we settled on pixel a long time ago. For the data processors you would think the user should be aware of what they want and select the registration they need. We may even consider adding a warning to any grid processor given a remove file without registration that they did not specify registration and we have defaulted to gridline registration.

Happy to hear comments on this before moving on it.

[PaulWessel]

Specifically, we have this in the gmt_init_module function:

if (!strcmp (mod_name, "grdtrack")) API->use_gridline_registration = true; /* Override API default since grdtrack is a data processor */

[joa-quim]

I have nothing against processors use grid reg, but at least for nc grid when they were not produced by us we guess registration from x,y coordinates and from GDAL we apply the rule if not explicit in file grids are grid and images are pix registered. This together with the native binary formats pretty much cover everything so we should have nothing with a missing registration info.

[PaulWessel]

Well, here we know exactly that things are since we produced these grids. But still need to decide, no?