Potential Region-Dependent Artifact: Vertical Striping Appears in RTP Results Under Specific Geomagnetic Parameters

[lazyn1997]

Problem Description

When using the harmonica to perform reduction to the pole (RTP) on magnetic data, noticeable vertical stripes (or smoothing artifacts) appear in the results. The occurrence of this artifact shows a strong correlation with the geomagnetic parameters (particularly magnetic inclination and declination) of the survey area.

Steps to Reproduce

1. Prepare two gridded magnetic datasets from different regions.
2. Region 1 (Problematic Case): Perform RTP using the following parameters:
   • Central coordinates: -66.2, -23.0
   • Magnetic declination: -4.3°
   • Magnetic inclination: -20°
3. Region 2 (Comparison Case): Perform RTP using the following parameters:
   • Magnetic declination: 6.4°
   • Magnetic inclination: -40.3°
4. Compare the RTP result maps for the two regions.
5. Observation: The RTP result for Region 1 (-20° inclination) exhibits clear vertical striping/smoothing artifacts. This phenomenon is almost invisible in the result for Region 2 (-40.3° inclination).

Expected Behavior

The RTP processing should yield results free of obvious directional artifacts regardless of the geomagnetic parameters. Both regions should exhibit consistent data quality.

    def _reduce_to_pole_harmonica(self, data, grid_lon, grid_lat, inclination, declination, improve=False):
        """
        使用Harmonica进行化极处理 - 修正版
        """
        # 1. 验证坐标顺序和单调性
        x_coords = grid_lon[0, :]  # 经度坐标（东向）
        y_coords = grid_lat[:, 0]  # 纬度坐标（北向）

        # print(f"坐标范围: X({x_coords.min():.2f} - {x_coords.max():.2f}), "
        #       f"Y({y_coords.min():.2f} - {y_coords.max():.2f})")
        flipud = False
        fliplr = False
        # 检查坐标是否单调递增
        if not (np.all(np.diff(x_coords) > 0) and np.all(np.diff(y_coords) > 0)):
            if not np.all(np.diff(x_coords) > 0):
                print("警告: X坐标不是严格单调递增，可能需要调整")
                # 如果坐标递减，需要反转
                if np.all(np.diff(x_coords) < 0):
                    x_coords = x_coords[::-1]
                    data = data[:, ::-1]
                    fliplr = True
                    print("X坐标已反转")
                else:
                    raise ValueError("X坐标不是严格单调递增且不是严格单调递减，请检查数据")
            if not np.all(np.diff(y_coords) > 0):
                print("警告: Y坐标不是严格单调递增，可能需要调整")
                if np.all(np.diff(y_coords) < 0):
                    y_coords = y_coords[::-1]
                    data = data[::-1, :]
                    flipud = True
                    print("Y坐标已反转")
                else:
                    raise ValueError("Y坐标不是严格单调递减且不是严格单调递增，请检查数据")
        else:
            print("X坐标已验证为严格单调递增且Y坐标已验证为严格单调递增")

        # 2. 创建正确的DataArray
        # 关键修正：确保维度顺序正确
        grid = xr.DataArray(
            data=data,
            dims=("northing", "easting"),  # 使用更明确的维度名称
            coords={
                "northing": y_coords,  # 北向坐标（纬度）
                "easting": x_coords,  # 东向坐标（经度）
            },
            attrs={
                "inclination": inclination,
                "declination": declination
            }
        )

        # inclination 地磁场的倾角
        # declination 地磁场的偏角
        # magnetization_inclination 异常体的总磁化倾角，设为 None，harmonica 会默认等于 inclination
        # magnetization_declination 异常体的总磁化偏角，设为 None，harmonica 会默认等于 declination
        rtp_data = hm.reduction_to_pole(
            grid,
            inclination=inclination,
            declination=declination,
            magnetization_inclination=inclination,
            magnetization_declination=declination
        )
        rtp_values = rtp_data.values
        if flipud:
            rtp_values = np.flipud(rtp_values)
            print("对化极结果进行上下翻转")
        if fliplr:
            rtp_values = np.fliplr(rtp_values)
            print("对化极结果进行左右翻转")

        def stabilize_fft(rtp, level=4):
            fft_rtp = fft.fft2(rtp)
            ny, nx = rtp.shape

            kx = np.fft.fftfreq(nx)
            ky = np.fft.fftfreq(ny)
            kx, ky = np.meshgrid(kx, ky)
            k = np.sqrt(kx ** 2 + ky ** 2)

            LP = np.exp(-(k ** 2) * level)  # Gaussian low-pass

            return np.real(fft.ifft2(fft_rtp * LP))

        def rtp_regularize(rtp, alpha=0.15):
            return rtp / (1 + alpha)

        def edge_taper(rtp, w=0.1):
            ny, nx = rtp.shape
            wx = np.hanning(int(nx * w))
            wy = np.hanning(int(ny * w))

            wx = np.concatenate([wx[:len(wx) // 2], np.ones(nx - len(wx)), wx[len(wx) // 2:]])
            wy = np.concatenate([wy[:len(wy) // 2], np.ones(ny - len(wy)), wy[len(wy) // 2:]])

            taper = wy[:, None] * wx[None, :]
            return rtp * taper

        def improve_harmonica_rtp(rtp):
            rtp1 = stabilize_fft(rtp, level=6)  # 关键滤波稳定方向
            rtp2 = rtp_regularize(rtp1, alpha=0.12)  # 抑制异常区域爆炸
            rtp3 = edge_taper(rtp2, w=0.08)  # 去除边界振铃
            return rtp3

        if improve:
            rtp_values = improve_harmonica_rtp(rtp_values)

        return rtp_values

harmonica result.docx

[santisoler]

Hi @lazyn1997. Thanks for opening this issue.

I understand that for the region that has an inclination of -20 degrees, the RTP shows noticeable stripes along the latitudinal direction, which are not desired. Nonetheless, the reason why those stripes appear can have multiple sources. For example, they might some artifact present in the grid created by aliasing, if no proper anti-aliasing technique was used when creating the grid. Also, an inclination of -20 degrees is expected to be somewhat problematic for computing RTP (since the inducing field is almost horizontal with respect to the Earth surface), so comparing it with the -40 degrees of inclination might not be a 100% one-to-one comparison.

That being said, I think using some kind of filter (like the low-pass and the Hanning filter) might help in removing those artifacts. For example, if the grid contains some artifacts due to aliasing, those filters might be removing them.

The RPT function we have in Harmonica implements the simplest definition of the reduction-to-the-pole method. I've taken a second look at it and the code seems correct. Not using any extra filter in the RTP function is a design choice. The good thing is that it can always be extended in case you want to use any other filter on top of that, like you did. You can even use the hm.filters.reduction_to_the_pole_kernel directly too.

Please, let me know what do you think!

[santisoler]

I've just opened #630 to extend the tests we have on that function, in case there's anything wrong in it that we might not be noticing.

[lazyn1997]

Thank you very much for your detailed explanation and for taking a second look at the RTP implementation in Harmonica. I really appreciate your insights.

In our tests, we performed RTP on the same study area using both the real geomagnetic inclination (−20°) and a simulated larger inclination. We found that increasing the inclination significantly reduces the latitudinal striping observed in the RTP results. This comparison is documented in compare.docx. This suggests that at least part of the striping is related to the low-inclination.

We also reviewed the Harmonica documentation and related references, where it is stated that, without additional treatment, RTP at low latitudes tends to amplify high-frequency noise.

Our understanding is that this amplification may correspond to the vertical striping patterns seen in the results. We experimented with several post-processing filters (e.g., ndimage.gaussian_filter and ndimage.median_filter), but the improvements were limited. In this context, we would like to ask whether Harmonica provides, or plans to provide, any recommended or more optimal post-processing strategies for RTP results at low latitudes.

In addition, we are aware that reduction to the pole based on the equivalent-source method is often considered a more stable solution for low-latitude magnetic data. We noticed from the repository discussions that an EquivalentMagSources approach is planned for v0.8.0 #426. Could you please share the current status of this development? Also, in the meantime, is it feasible to approximate this method in the current release by combining existing Harmonica functionalities (e.g., equivalent sources + forward modeling)?

Finally, from your perspective, which approach would you consider more appropriate for addressing this issue in practice:
(1) applying RTP with carefully designed filtering, or
(2) using an equivalent-source-based RTP approach?
If possible, could you provide some guidance or example code? I would also be happy to share the original GRD file as well as the Oasis Montaj RTP result if that would be helpful.

Thanks again for your time and for maintaining such a useful library.

compare harmonica result.docx