Try to get matrix as grid to work with changing central meridian

[PaulWessel]

See #3811 for background. This PR extends the test to give three global plots for different central longitudes: 0, 98.7, and 180. I am adding the correct PS original for this test as well.

However, this effort runs into trouble with the function gmt_grd_setregion that determines what region of the grid we should use. As you can see in the Changes fo that function, we are looking at the "easy" case of a global grid. Turns out not to be easy. We compute a w/e extent and then check if it exceeds or equals 360. If it does we trim back one gridnode at the time. In the case of the matrix feed, we end up with a 359 range instead and plot fails to rotate the longitude (see lines 5485 in gmt_api.c in Changes). All my efforts in trying to have it become 360 causes other tests to fail, such as globalgrids.sh. Right now, the apimat_360.sh test fails for the longitude of 98.7. It will work fine for any central meridian that is a multiple of the increment, here 1 degree. If I tinker with the true_global_region parameter by removing the check for registration then apimat_360.sh works, but then globalgrid.sh fails, and vice versa.

I am leaving this as WIP for now. I will try to run the equivalent 98.7 longitude with the grid (which works fine) and see if I can learn what the hell is different here.

[seisman]

@PaulWessel Do you want to approve this PR and merge it into master?

[PaulWessel]

No, you can wait until things work. Running

gmt grdimage @earth_relief_01d -JH98.75/6i -B -pdf map

I find that the gmt_grd_setregion determines the region to be -81/278, in other words 359 degrees. However, since there is no check on rotating here (since we read it in that way from file) the only thing that happens is the projection to inches for plotting uses 180 rows (correct) and 359 instead of 360 columns. Plot looks OK since you cannot tell 360 from 359 columns, but it is wrong. Yet, changing that part of the code (above) to find 360 messes up tests.

[PaulWessel]

I think I need to apply the fix so that this is 360 and then see what goes wrong with the globalgrid.sh script.

[PaulWessel]

Hi @seisman, I think I have reached a happy solution on this. By pursuing the possibility that the globalgrid.sh test might not be correct given the evidence that the changes discussed above actually resolved the correct domain and dimensions of the global relief grid, I found that the globalgrid.sh test indeed had influenced the coding too much. Of the 8 panels in that test, only the lower left "fails" with the proposed changes. Before my changes it passed and looked like this:

This is a global, pixel grid defined on a 0/360/-90/90 lattice with 30 degree spacing (12x6 nodes). It is plotted on a map using -R-165/195/-90/90 whose boundary cuts through a column of pixels. The old system truncates the right side of the grid at 180 instead of going to 210 outside the right boundary and showing a half pixel column. My proposed changes in gmt_grdio.c caused a problem in that the east boundary remained forced at 180 but it now counted all 12 pixels, thus fitting one extra node between 150W and 180 so that each pixel extended less than 30 degrees. To compensate, I made the proposed adjustment to grdimage where in a local function it computes the "projected" region in the case of no map projection. Since in this case we have a geographic grid with periodic longitudes, the last pixel x coordinate wraps around to the beginning and hence it was not used to indicate the full width of that image. By checking for this special case (wrap around in a Cartesian plot) I now show all 12 nodes with one clipped at the east boundary. Thus the new plot looks like this:

which is 50% better I think. Now, you may say why can't we also paint the beginning from 165W to 150W as well since it is a global grid. Part of the reason is that such a process would involve adding an extra column to the grid by duplicating the right column so it goes from 180W to 210E for a 390 degree range. That introduces many complications. For most data sets this is not going to be a problem since they tend to be larger than 12x6 nodes. I agree though that GMT should aim for perfection but we have never dealt with adding a column because of fractional coverage. So, since this PR addresses the original problem reported in #3811 I think this half-victory for globalgrid.sh is a bonus. I can make a feature request to deal with this problem in the future.

[seisman]

@PaulWessel The PS file test/api/apimat_360.ps is wrong.

map-matrix.png is generated by passing matrix to GMT (i.e., by converting apimat_360.ps to PNG), map-grid.png is generated by running equivalent GMT commands.

Comparing these two figures, you will see that in the top two panels, the west hemisphere are shifted a little.

[PaulWessel]

Yes, you are right. The wrapping is not quite right. Sigh. Dinner time.

[PaulWessel]

Trying one more time, this with an obvious fix to the longitude issue.

[seisman]

The test looks correct now.

The PyGMT case (GenericMappingTools/pygmt#515) works well using GMT_IN or GMT_IN|GMT_IS_DUPLICATE, but is still wrong (sometimes crashes) using GMT_IN|GMT_IS_REFERENCE.

[PaulWessel]

OK, my testapi_matrix_360 is passing the matrix as double so of course the REFERENCE is overwritten to be DUPLICATE, so I am not testing the REFERENCE case. So I will change to float so we can test REFERENCE since DUPLICATE works

[PaulWessel]

Since this works for duplciate, lets approve this so I can start another that adds the proper reference. I need to make a few other changes first, @seisman. The only difference since a few hours ago is line 5495 in gmt_api.c that restricted REFERENCE to external memory. I do not think that is correct in general.