Fix north-south flip in FFT filters for descending-coordinate grids

src/harmonica/filters/_utils.py

@@ -79,6 +79,18 @@ def apply_filter(
     # xrft doesn't know what to do with them.
     non_dim_coords = {c: grid[c] for c in grid.coords if c not in grid.indexes}
     grid = grid.drop_vars(non_dim_coords.keys())
+    # The FFT-based filters assume the grid coordinates are in ascending order:
+    # xrft builds the frequency coordinates from positive sample spacings, and
+    # the inverse transform always returns the grid in ascending order. Grids
+    # read from some file formats store the northing in descending order, which
+    # silently produced a north-south flipped output because the original
+    # (descending) coordinates were assigned to the ascending result (see #586).
+    # Sort the grid to ascending order before transforming and restore the
+    # original coordinate order at the end.
+    original_coords = {dim: grid[dim] for dim in dims}
+    needs_reorder = any(bool(grid[dim][0] > grid[dim][-1]) for dim in dims)
+    if needs_reorder:
+        grid = grid.sortby(list(dims))
     if pad:
         # By default, use a padding width of 25% of each grid dimension.
         # Fedi et al. (2012; doi:10.1111/j.1365-246X.2011.05259.x) suggest
@@ -109,6 +121,12 @@ def apply_filter(
     filtered_grid = filtered_grid.assign_coords(
         {dims[1]: grid[dims[1]], dims[0]: grid[dims[0]]}
     )
+    # Restore the original coordinate order (e.g. descending northing) so the
+    # output is aligned with the input grid and not flipped (see #586).
+    if needs_reorder:
+        filtered_grid = filtered_grid.reindex(
+            {dim: original_coords[dim] for dim in dims}
+        )
     # Restore the non-dimensional coordinates if desired
     if not drop_coords:
         filtered_grid = filtered_grid.assign_coords(

test/test_transformations.py

@@ -301,6 +301,32 @@ def test_derivative_upward(sample_potential, sample_g_z):
     assert rms / np.abs(g_up).max() < 0.015
 
 
+@pytest.mark.parametrize("flip_dim", ["northing", "easting"])
+def test_filter_descending_coordinate_not_flipped(sample_potential, flip_dim):
+    """
+    Filters should not flip the output for grids with descending coordinates.
+
+    Grids read from some file formats store the northing (or easting) in
+    descending order. The FFT-based filters used to assign the original
+    descending coordinates to the ascending result of the inverse transform,
+    silently flipping the output along that dimension (see
+    https://github.com/fatiando/harmonica/issues/586).
+    """
+    ascending = derivative_upward(sample_potential)
+    # Same physical grid, but with one coordinate stored in descending order
+    flipped_grid = sample_potential.isel({flip_dim: slice(None, None, -1)})
+    assert flipped_grid[flip_dim].values[0] > flipped_grid[flip_dim].values[-1]
+    descending = derivative_upward(flipped_grid)
+    # The output must keep the descending coordinate order of its input
+    assert descending[flip_dim].values[0] > descending[flip_dim].values[-1]
+    # Sorting both results on the flipped dimension must make them identical:
+    # if the output were flipped, the values would not match after sorting.
+    npt.assert_allclose(
+        ascending.sortby(flip_dim).values,
+        descending.sortby(flip_dim).values,
+    )
+
+
 def test_derivative_upward_order2(sample_potential, sample_g_zz):
     """
     Test higher order of derivative_upward function against the sample grid.