Allow for grids with negative lat/lon increments

[MarkWieczorek]

pygmt assumes that both the latitudes and longitudes of an xarray Dataset increase with increasing indices of the array. However, it is common for the latitudes to be arranged with decreasing values such that the first row corresponds to 90 N and the last to 90 S. Furthermore, in some weird and very old cases, planetary longitude has been defined as increasing in the westward direction.

Here, I have modified the function dataarray_to_matrix to check whether the longitude and latitude increments are positive. If they are negative, the array is flipped on the corresponding dimension.

There are four possible scenarios, and I have attempted to avoid creating an additional array. If there is a better way to do this, please let me know!

Fixes # #368

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to docstrings or tutorials.

[welcome]

💖 Thanks for opening this pull request! 💖

Please make sure you read our contributing guidelines and abide by our code of conduct.

A few things to keep in mind:

• If you need help writing tests, take a look at the existing ones for inspiration. If you don't know where to start, let us know and we'll walk you through it.
• All new features should be documented. It helps to write the docstrings for your functions/classes before writing the code. This will help you think about your code design and results in better code.
• No matter what, we are really grateful that you put in the effort to do this! 🎉

[weiji14]

Right, I've taken a hard look at this, and you've definitely done a good job at getting the functionality right! I think we can improve the readability of the code though, e.g. by using np.flip instead of grid.values[:, ::-1]. Edit: Found out that xarray has a sortby function which means we can simplify all of the if-then statements into one line. Will recommend some suggested changes shortly.

I'm not sure if you're familiar with writing tests, but could you add some for the dataarray_to_matrix function? Just add them above these ones:

pygmt/pygmt/tests/test_clib.py

Lines 767 to 786 in 55af9b6

	def test_dataarray_to_matrix_dims_fails():
	"Check that it fails for > 2 dims"
	# Make a 3D regular grid
	data = np.ones((10, 12, 11), dtype="float32")
	x = np.arange(11)
	y = np.arange(12)
	z = np.arange(10)
	grid = xr.DataArray(data, coords=[("z", z), ("y", y), ("x", x)])
	with pytest.raises(GMTInvalidInput):
	dataarray_to_matrix(grid)
	def test_dataarray_to_matrix_inc_fails():
	"Check that it fails for variable increments"
	data = np.ones((4, 5), dtype="float64")
	x = np.linspace(0, 1, 5)
	y = np.logspace(2, 3, 4)
	grid = xr.DataArray(data, coords=[("y", y), ("x", x)])
	with pytest.raises(GMTInvalidInput):
	dataarray_to_matrix(grid)

To get you started on the unit tests, I recommend constructing a diagonal array like this, which you can flip around in the x/y directions by changing the start/stop values:

    data = np.diag(v=np.arange(3))
    x = np.linspace(start=4, stop=0, num=3)
    y = np.linspace(start=5, stop=9, num=3)
    grid = xr.DataArray(data, coords=[("y", y), ("x", x)])

Just give me a shout out if you need any help 😄

[weiji14]

So this uses xarray's builtin sort mechanism that simplifies the code by heaps.

Suggested change

	if inc[1] < 0:
	if any([i < 0 for i in inc]): # Sort grid properly if there are negative increments
	inc = [abs(i) for i in inc]
	grid = grid.sortby(variables=list(grid.dims), ascending=True)
	matrix = as_c_contiguous(grid.values[::-1])

Note ⚠️: do not commit this change directly, you'll need to manually edit the file and make the change to the lines accordingly.

Also, it would be good to see some unit tests to make sure it actually works for every possible combination of flipped latitude/longitude coordinates.

[MarkWieczorek]

That works. I suspect that the sorting operation is slower than simply flipping the array, but compared to everything involved with the projection routines, this is certainly not a limiting operation.

As for the tests, I'd be happy to try, but I'm not sure what the point would be. It seems that the test would be to check if xarray works properly, not pygmt. What kind of test were you thinking of?

Notes

• It appears that this works even if the DataSet doesn't supply coordinates (the input grid is unmodified). Though this probably won't occur in pygmt.
• The only problem I could envision is if the user defined the longitude coordinates as [0, 90, 180, 270, 0]. In this case the sort operation would provide incorrect results. It is hard to imagine someone doing this. If there is a test for this scenario, I would probably redefine the last coordinate somewhere else in the code. As part of such a test, you might also want to check that the two data values are the same (I think that gmt might already do this). This scenario is caught when testing if the coordinate increments are all the same or not.

[weiji14]

If speed is the concern, it would be best to pass in a file anyway. You're right that the test is more of a check on xarray than pygmt, but we're actually not testing that dataarray_to_matrix properly now so it's a good time as any to add it.

[leouieda]

If speed is the concern, it would be best to pass in a file anyway.

I'm not sure even GMT handles this properly. So passing in a file probably doesn't help much.

[weiji14]

@leouieda, do you want to take a look at this dataarray_to_matrix function if you have time?

[leouieda]

As for the tests, I'd be happy to try, but I'm not sure what the point would be. It seems that the test would be to check if xarray works properly, not pygmt. What kind of test were you thinking of?

@MarkWieczorek for the test, the thinking to check that if we pass in an array with flipped latitude or longitude that dataarray_to_matrix really does flip the matrix as expected (like what @weiji14 suggested). So we're not checking xarray, but checking that we wrote the sortby correctly (and that future changes don't break the functionality).

[weiji14]

Just checking in here @MarkWieczorek. Have you had a chance to write the tests yet? I'd be keen to merge this PR as soon as possible as I've been hitting the same flipped lat/lon problem a lot recently too 😂

[MarkWieczorek]

Sorry but no. I all likelihood, this will probably take a week or two (I'm busy, plus I've never written pytests before). If you do this for me, I promise instead to improve the web documentation later...

Also, I think that this issue is probably related: #375

I'm not sure exactly what is going on, but I think that dataarray_to_matrix probably needs to check if the last redundant longitude column is missing, and if so, to add it.

[weiji14]

Deal, I'll write up the tests for this now. I did take a look at that other issue, but it will require a bit more investigation.

[welcome]

🎉🎉🎉 Congrats on merging your first pull request and welcome to the team! 🎉🎉🎉

Please open a new pull request to add yourself to the AUTHORS.md file. We hope that this was a good experience for you. Let us know if there is any way that the contributing process could be improved.

[weiji14]

Thanks @MarkWieczorek! Hopefully this resolves your original issue at #368, I know it'll definitely help my workshop tutorials run more smoothly tomorrow. Definitely feel free to contribute more in the future, I can see you've got a lot of great ideas 😄

[leouieda]

Thank you for the contribution Mark! 🎉