Missing parameters when partitioning spectra

[DrakonianMight]

Hi,
I am trying to understand the behaviour for the generation of parameters when I am partitioning the spectrum.
So my case is that I am wanting to partition the spectrum down into quite a narrow frequency range (which I know can cause issues, with reintegration), but I didn't expect there to be no peak parameters produced.

Here's my code:

myspec = wavespectra.read_ww3("myfile")
st = 12
ed = 19
part = myspec.spec.split(fmin=1/ed, fmax =1/st ).chunk({"freq": -1})
part.spec.tp(smooth = False).plot()

And here is a plot of the tp.

I have had a quick look at the source code, but there appears to be a couple of different functions that calculate the peak parameters and I am not sure which one is currently in use. My aim here is to have a tp value produced, even if it at the boundary of the frequency range that I have supplied (the same with other peak parameters like peak direction).

Cheers

[rafa-guedes]

@DrakonianMight the peak-based parameters such as dpm or tp require a local peak in the 1d frequency spectrum. If no local peak exists those functions return nans. It is not uncommon to end up with partitions that have no peak when we split the spectrum based on an arbitrary frequency threshold and I assume this is the cause of your missing values. You can check that by plotting the 1d spectra for the partitions for which you have identified missing values. Notice that some stats that do not require the frequency location of the peak to be defined (such as hs()) still return values.

[DrakonianMight]

Thanks @rafa-guedes, yeah I can see that for the partitions in question. Yeah, I can see that. The behaviour I am looking for is to be able to take the peak at the boundary, if a peak is not present in the partitioned spectrum.

I wondered whether there was a way of parsing the ipeak in based on just using some sort of argmax.

[rafa-guedes]

In the beginning we used argmax which is faster and a lot simpler to calculate since it can be done directly with dask and does not require those ufuncs we needed to write to allow indexing using ipeak. However in cases like this example it isn't really capturing a peak hence we moved to the new approach that requires local maxima which we thought was more physically sound (one could argue though that an arbitrary frequency split is also not physically sound). But unless there is a strong argument to do it I'd prefer to keep it the way it is now. It should be easy enough though to adjust your code so you can use the maximum frequency of the low-frequency partition (or the minimum of the high-frequency one) by doing something like:

tp = part.spec.tp(smooth=False).fillna(1 / part.freq.max())

I haven't tested this snippet but something along these lines should work.

[DrakonianMight]

Yeah I absolutely agree with the thinking. In my experience the use case comes up a fair amount about splitting the spectrum using more than just the topographical/watershedding method. Like you say, I am not sure how "valid" the spectrum is, particularly when the number of frequency bins in the model may be quite limited. Perhaps consider it as an optional call?

Thanks for the snippet for Tp, thats a neat one liner (I was writing something far more complex). I guess the problem still remains about calculating peak direction, any tips?

[rafa-guedes]

True that will only solve it for tp.

Just so you understand it, the spectral stats requiring the peak frequency index ipeak are defined as ufuncs in wavespectra.core.xrstats so they can be calculated lazily with dask. Each ufunc uses numpy-based core functions (written with numba for speed) that are defined in wavespectra.core.npstats. You could check out the core functions you need (such as dpm) and write some custom functions to calculate it using ipeak as the index of the boundary.

If you really want this as a functionality in wavespectra you could try to implement it as an optional argument inside each ipeak-based method, something like fill_value="max" which would take NaN by default but could use argmax if max was chosen. I'd be open to a pull request on this if you want to have a go.