Bug in matlab function shapiro1

[bjorn]

Hi,

It seems to be a bug in the matlab script utility/shapiro1.m

Higher order filtering does not work properly.
The supposedly most used 2nd order is OK as far as I can see.

The loop overwrites cor(2:Im-1) each time,
making cor(2:Im-1)=2.0.*Finp(2:Im-1) - Finp(1:Im-2) - Finp(3:Im)
independent of the order.
In particular, the stencil is always three points wide
(it should be five for 4th order, ...)

The bug is also present in the python translation in
Kate Hedstrom's pyroms package on github.

Bjørn Ådlandsvik

[arango]

Thank you. I will look at this problem when I get the chance. If you corrected this function, please send me the script so I can take at look at it.

[bjorn]

Hi,

I came to this routine by tracing backwards from the python version
in pyroms. Unfortunately, I do not have a working high order version.
And, perhaps due to the bug, I have never really understood what the
function is trying to do.

The point in Shapiro 1970-paper is to make higher order by successive
filtering with simple 3-point stencils. So at least the loop should
provide successive corrections. Neglecting the boundary second order
is computed by

cor(2:Im-1) = 2.0*Finp(2:Im-1) - Finp(1:Im-2) - Finp(3:Im);
F2 = Finp - 0.25*cor;

The next step can then be done by

cor(2:Im-1) = 2.0*F2(2:Im-1) - F2(1:Im-2) - F2(3:Im);
F4 = F2 + 0.25*cor;

The stencils are respectively [1 2 1]/4 and [-1 4 10 4 -1]/16
as described in the Shapiro papers. After this it gets more
complicated.

For second order the matlab function is correct at open sea, but
without mask treatment the results may be polluted by land values.
The shapiro filter in the ROMS code (Utility/shapiro.F) is only
second order but has improved land handling by doing nothing in sea
cells adjacent to land cells.

I suppose the higher order version is not used very much. That may be
the reason the bug has not been noticed. One possibility is to stay
consistent with the ROMS code and provide only second order but with
proper land treatment. I am not sure if the ROMS' way is the best
way to treat land.

Bjørn