From 6bfdd66f643cadd77ba59ec467556f83258849a7 Mon Sep 17 00:00:00 2001
From: gaoflow <gaobing1230@gmail.com>
Date: Mon, 1 Jun 2026 18:13:23 +0200
Subject: [PATCH] Fix north-south flip in FFT filters for descending-coordinate
 grids

The FFT-based filters (derivatives, upward continuation, reduction to pole,
etc.) silently produced an output flipped along a dimension when the input
grid stored that coordinate in descending order, as happens for grids read
from some file formats (e.g. ERMapper/ERS) where northing decreases with row
index.

The inverse FFT always returns the grid in ascending coordinate order, but
apply_filter then assigned the original (descending) coordinates onto that
ascending result, mislabelling every row/column and flipping the data
relative to its coordinates. The frequency coordinates built by xrft also
assume positive sample spacing.

Sort the grid to ascending order before transforming and restore the original
coordinate order afterwards via reindex. Grids already in ascending order are
unaffected. Adds a parametrised regression test (northing and easting)
asserting a descending-coordinate grid is not flipped.

Fixes #586
---
 src/harmonica/filters/_utils.py | 18 ++++++++++++++++++
 test/test_transformations.py    | 26 ++++++++++++++++++++++++++
 2 files changed, 44 insertions(+)

diff --git a/src/harmonica/filters/_utils.py b/src/harmonica/filters/_utils.py
index 5971e4558..de7dc81c3 100644
--- a/src/harmonica/filters/_utils.py
+++ b/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(
diff --git a/test/test_transformations.py b/test/test_transformations.py
index 19209f6f3..20f9d7eb2 100644
--- a/test/test_transformations.py
+++ b/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.