Add special check for non-rotated global grids

[PaulWessel]

We usually read global grids from file and they get rotated depending on the wesn selection. However, when a global grid is passed from another environment, as GMT_IS_REFERENCE, we may be projecting that image in gmt_grd_project and if the repeated meridians fall inside the map (since no rotation has happened) they are counted twice, yielding a slightly different result. Only applies to gridline-registrered grids. This fix detects this special situation and skips one of the repeated meridians. Addresses the initial issue raised in #3844.
Note sure about backport. Was the gmt_whole_earth backported?

[seisman]

This PR should be backported.

Here are diff images generated by matplotlib image comparison functions. I find then more readable than the gm diff images.

This is from GMT master branch:

This is from this PR:

The central meridian is fixed by this PR, but please see if you can reduce the right side differences.

[PaulWessel]

Not sure. It was relatively clear what to do for when the repeated meridian lands inside. While I am not sure why there is nothing on the west border, I would suspect this has to do with BCs. When the grid is read in and we know it is global, we use actual data values from the other side to set the two extra rows/cols in the grid data. THe matrix has no padding, so when it is passed into gmt_grd_project, this happens:

/* Only input grid MUST have at least 2 rows/cols padding - otherwise we must allocate a temp grid */
if (I->header->pad[XLO] < 2 || I->header->pad[XHI] < 2 || I->header->pad[YLO] < 2 || I->header->pad[YHI] < 2) {
	unsigned int pad2[4] = {2, 2, 2, 2};
	GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "gmt_grd_project: Input grid has insufficient padding - create and work on a duplicate with r row/col pad\n");
	if ((I2 = GMT_Duplicate_Data (GMT->parent, GMT_IS_GRID, GMT_DUPLICATE_DATA, I)) == NULL) {
		GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmt_grd_project: Unable to duplicate grid\n");
		return GMT_RUNTIME_ERROR;
	}	
	gmt_grd_pad_on (GMT, I2, pad2);	/* Add pad */
	gmt_BC_init (GMT, I2->header);	/* Initialize grid interpolation and boundary condition parameters */
	gmt_grd_BC_set (GMT, I2, GMT_IN);	/* Set boundary conditions */
	I = I2;	/* Use this input grid instead */
}

So we end up duplicating the non-padded matrix so we can project. So assuming the grid/matrix is flagged as geographic, it should end up applying periodic boundary conditions at west and east and it should be the same for the grid as well as the matrix. But clearly the east boundary differs. Will see if I can find any inconsistencies.

[PaulWessel]

Think I have the answer: It is the opposite of the fix I added. For this projection and a gridline-registered grid, the longitude that corresponds to the map west and east will be plotted twice. Since the rotated grid (read from file) has a min/max that equal east/west, it is plotted twice and all is well.
For the GMT_IS_REFERENCE grid that is not rotated, the repeated -180 and +180 meridian both project to the the interior of the image (unless we use -Rd). So the earlier fix was to only project one of them and that fixed the reported problem. However, in the above case, w=-60 and e=300 but the duplicated longitude is still -180/180 for a referenced grid. That means that when we plot the single meridional line for west=-60 we must actually add 360 to its longitude and do it again - otherwise there is nothing projected at the east boundary. Will try to add such a fix.

[PaulWessel]

So I would like to debug my new algorithm but now I cannot get the python script to pass control to Xcode when it calls fig.grdimage(grid, projection="W120/15c", cmap="geo"). Despite making sure I have GMT_LIBRARY_PATH=/Users/pwessel/GMTdev/gmt-dev/xbuild/src/Debug as well as GMT_SHAREDIR=/Users/pwessel/GMTdev/gmt-dev/dbuild/gmt6/share, I run conda activate pygmt, then open up the python terminal. attach Xcode to the process (Python, the only one) and set stop points in Xcode in gmt.c, it just finishes. I have rebooted, rebuilt, etc. No go. This worked yesterday, so I must be missing something. My steps after reboot:

1. export GMT_LIBRARY_PATH=/Users/pwessel/GMTdev/gmt-dev/xbuild/src/Debug
2. export GMT_SHAREDIR=/Users/pwessel/GMTdev/gmt-dev/dbuild/gmt6/share
3. conda activate pygmt
4. python
5. Run the first few commands (import pygmt, etc)
6. Attach Xcode to process python.
7. Run fig.grdimage(grid, projection="W120/15c", cmap="geo")

but it just completes. Plot is made, but I am unable to debug. Any ideas?

[seisman]

set stop points in Xcode in gmt.c

Do you mean "gmt_api.c"?

[PaulWessel]

Hm, perhaps that is the problem. I am still in my CLI module being called by gmt but I need to think GMT_Call_Module and put a stop point inside Call_Module.

[PaulWessel]

All is well. Another lesson learned on debugging. 99% of my debugging has been via gmt calling GMT_Call_Modules so that has been my entry point into the library...

[PaulWessel]

Please give this another shot. My gm compare show zero difference between the two gridline plots and all tests pass (well minus the gdalnn and the imagepluscpt).

[seisman]

Great!