Allow for grids with negative lat/lon increments (#369)

Flips grids with negative latitude/longitude increments properly using xarray sortby(). Includes tests that check to make sure that dataarray_to_matrix works for standard case where longitude (x) and latitude (y) dimensions of grids are in positive increments (ascending order). Also ensures that grids are flipped correctly when the x and/or y axis are in negative increments (descending order).

pygmt/clib/conversion.py

@@ -9,13 +9,15 @@
 
 def dataarray_to_matrix(grid):
     """
-    Transform a xarray.DataArray into a data 2D array and metadata.
+    Transform an xarray.DataArray into a data 2D array and metadata.
 
     Use this to extract the underlying numpy array of data and the region and
     increment for the grid.
 
     Only allows grids with two dimensions and constant grid spacing (GMT
-    doesn't allow variable grid spacing).
+    doesn't allow variable grid spacing). If the latitude and/or longitude
+    increments of the input grid are negative, the output matrix will be
+    sorted by the DataArray coordinates to yield positive increments.
 
     If the underlying data array is not C contiguous, for example if it's a
     slice of a larger grid, a copy will need to be generated.
@@ -102,6 +104,11 @@ def dataarray_to_matrix(grid):
             )
         region.extend([coord.min(), coord.max()])
         inc.append(coord_inc)
+
+    if any([i < 0 for i in inc]):  # Sort grid when 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])
     return matrix, region, inc
 

pygmt/tests/test_clib.py

@@ -764,6 +764,58 @@ def test_write_data_fails():
                 )
 
 
+def test_dataarray_to_matrix_works():
+    "Check that dataarray_to_matrix returns correct output"
+    data = np.diag(v=np.arange(3))
+    x = np.linspace(start=0, stop=4, num=3)
+    y = np.linspace(start=5, stop=9, num=3)
+    grid = xr.DataArray(data, coords=[("y", y), ("x", x)])
+
+    matrix, region, inc = dataarray_to_matrix(grid)
+    npt.assert_allclose(actual=matrix, desired=np.flipud(data))
+    npt.assert_allclose(actual=region, desired=[x.min(), x.max(), y.min(), y.max()])
+    npt.assert_allclose(actual=inc, desired=[x[1] - x[0], y[1] - y[0]])
+
+
+def test_dataarray_to_matrix_negative_x_increment():
+    "Check that dataarray_to_matrix returns correct output with flipped x dimensions"
+    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)])
+
+    matrix, region, inc = dataarray_to_matrix(grid)
+    npt.assert_allclose(actual=matrix, desired=np.flip(data, axis=(0, 1)))
+    npt.assert_allclose(actual=region, desired=[x.min(), x.max(), y.min(), y.max()])
+    npt.assert_allclose(actual=inc, desired=[abs(x[1] - x[0]), abs(y[1] - y[0])])
+
+
+def test_dataarray_to_matrix_negative_y_increment():
+    "Check that dataarray_to_matrix returns correct output with flipped y dimensions"
+    data = np.diag(v=np.arange(3))
+    x = np.linspace(start=0, stop=4, num=3)
+    y = np.linspace(start=9, stop=5, num=3)
+    grid = xr.DataArray(data, coords=[("y", y), ("x", x)])
+
+    matrix, region, inc = dataarray_to_matrix(grid)
+    npt.assert_allclose(actual=matrix, desired=data)
+    npt.assert_allclose(actual=region, desired=[x.min(), x.max(), y.min(), y.max()])
+    npt.assert_allclose(actual=inc, desired=[abs(x[1] - x[0]), abs(y[1] - y[0])])
+
+
+def test_dataarray_to_matrix_negative_x_and_y_increment():
+    "Check that dataarray_to_matrix returns correct output with flipped x/y dimensions"
+    data = np.diag(v=np.arange(3))
+    x = np.linspace(start=4, stop=0, num=3)
+    y = np.linspace(start=9, stop=5, num=3)
+    grid = xr.DataArray(data, coords=[("y", y), ("x", x)])
+
+    matrix, region, inc = dataarray_to_matrix(grid)
+    npt.assert_allclose(actual=matrix, desired=np.fliplr(data))
+    npt.assert_allclose(actual=region, desired=[x.min(), x.max(), y.min(), y.max()])
+    npt.assert_allclose(actual=inc, desired=[abs(x[1] - x[0]), abs(y[1] - y[0])])
+
+
 def test_dataarray_to_matrix_dims_fails():
     "Check that it fails for > 2 dims"
     # Make a 3D regular grid