search

MHKiT-Software / MHKiT-Python

MHKiT-Python provides the marine renewable energy (MRE) community tools for data processing, visualization, quality control, resource assessment, and device performance.
https://mhkit-software.github.io/MHKiT/
BSD 3-Clause "New" or "Revised" License
47 stars 45 forks source link

Contour Sampling - interpolation errors #218

Closed Graham-EGI closed 7 months ago

Graham-EGI commented 1 year ago

Describe the bug:

Sampling contours when using the full TP range to generate samples can throw errors about out of bounds.

To Reproduce:


from mhkit.wave.contours import samples_contour
import numpy as np

contour = {
        "hm_contour": [
            6.5104375921225985,
            6.481605176051669,
            6.395878699038057,
            6.25553964688066,
            6.064289128294204,
            5.827101053777997,
            5.550026500980775,
            5.2399588762209115,
            4.904370918721658,
            4.551036121492375,
            4.187747517321664,
            3.822046619562206,
            3.460974316011837,
            3.1108536945515706,
            2.777112323837948,
            2.464148735361109,
            2.1752455638505004,
            1.9125310930118724,
            1.6769925007134094,
            1.468546676049763,
            1.2861740867858051,
            1.1281132010528732,
            0.9920975562553997,
            0.8756026195093235,
            0.7760654957649292,
            0.6910506813219514,
            0.6183529334464551,
            0.5560441429858728,
            0.5024793166065582,
            0.4562776955682307,
            0.41629184455826873,
            0.38157326274604525,
            0.35133938235114825,
            0.3249442449464207,
            0.3018535876514179,
            0.28162424324031116,
            0.26388738502821346,
            0.24833502895234727,
            0.23470921487534,
            0.22279335404298675,
            0.21240531135348048,
            0.20339187113274412,
            0.1956243058058266,
            0.18899482605981446,
            0.18341373922113793,
            0.17880718104952067,
            0.17511531667298555,
            0.17229093056527775,
            0.17029834470479974,
            0.16911261952813073,
            0.16871900496696582,
            0.16911261952813073,
            0.17029834470479974,
            0.17229093056527792,
            0.17511531667298555,
            0.17880718104952067,
            0.18341373922113793,
            0.18899482605981446,
            0.1956243058058266,
            0.20339187113274412,
            0.21240531135348048,
            0.22279335404298675,
            0.23470921487534,
            0.24833502895234727,
            0.26388738502821346,
            0.28162424324031116,
            0.3018535876514179,
            0.3249442449464214,
            0.35133938235114825,
            0.38157326274604525,
            0.41629184455826873,
            0.45627769556823217,
            0.5024793166065582,
            0.5560441429858728,
            0.6183529334464551,
            0.6910506813219514,
            0.7760654957649292,
            0.8756026195093235,
            0.9920975562553997,
            1.1281132010528732,
            1.2861740867858051,
            1.468546676049763,
            1.6769925007134094,
            1.9125310930118724,
            2.1752455638505004,
            2.464148735361109,
            2.777112323837948,
            3.1108536945515706,
            3.460974316011837,
            3.822046619562206,
            4.187747517321664,
            4.551036121492375,
            4.904370918721658,
            5.2399588762209115,
            5.550026500980775,
            5.827101053777997,
            6.064289128294204,
            6.25553964688066,
            6.395878699038057,
            6.481605176051669
        ],
        "period_contour": [
            9.117702088877243,
            9.40062833834292,
            9.697189276579449,
            10.013193904712036,
            10.355297864300931,
            10.73119738040044,
            11.149858767818301,
            11.621778662694853,
            12.15926275331467,
            12.776697353689858,
            13.490766140867976,
            14.320530773293937,
            15.287247535020407,
            16.413736308907957,
            17.723068219565135,
            19.236329255304998,
            20.969307589216804,
            22.92820774668656,
            25.104932344799064,
            27.472974915885082,
            29.985249384370768,
            32.574928422818964,
            35.15953079263061,
            37.64747859097976,
            39.9457069013719,
            41.966920489565545,
            43.63556263212035,
            44.89210681625387,
            45.69565571740273,
            46.025003083118314,
            45.87836402616316,
            45.27198907071161,
            44.237887203209176,
            42.82089915579983,
            41.075375399938,
            39.061713852282736,
            36.84299477697488,
            34.48191497108843,
            32.0381743527164,
            29.566412014470544,
            27.11473252494289,
            24.72381268403931,
            22.426538275087083,
            20.248091824044977,
            18.20639631766962,
            16.312815099008954,
            14.57301261286454,
            12.98789170425094,
            11.554538134153567,
            10.267119512996755,
            9.117702088742424,
            8.096963444695891,
            7.194791381001411,
            6.400768752318263,
            5.704550842540303,
            5.096146242969578,
            4.566114563725966,
            4.105695099642138,
            3.706880230884472,
            3.3624462554180097,
            3.0659528469492345,
            2.8117206558050247,
            2.594794898857716,
            2.410901234715695,
            2.2563988590615316,
            2.128234610840843,
            2.0239009520855316,
            1.9413999476464496,
            1.879214777989282,
            1.8362897916733405,
            1.8120195247584943,
            1.8062462968045374,
            1.8192646559984538,
            1.8518286905022523,
            1.9051545657956632,
            1.9809052084298304,
            2.081137068997304,
            2.208182180879666,
            2.3644368825623765,
            2.5520391112056044,
            2.7724462222836355,
            3.025973402389324,
            3.3114007334623596,
            3.625773645224836,
            3.96448433151999,
            4.321640073688144,
            4.6906376678461905,
            5.064812168055523,
            5.438028735407125,
            5.805126408852355,
            6.162177189411112,
            6.506571227268509,
            6.836968084966102,
            7.153164226751702,
            7.455923264358824,
            7.746804800477568,
            8.028015463273823,
            8.302295169010733,
            8.572844048973165,
            8.843290936461374
        ]
}

# sample contour via number of samples across the line
min_tp = min(contour["period_contour"])
max_tp = max(contour["period_contour"])
# ten samples
period_samples = np.linspace(min_tp, max_tp, 10)
# error thrown about being above the interpolation range
hm_samples = samples_contour(period_samples, np.array(contour["period_contour"]), np.array(contour["hm_contour"]))

Minimal working example: A minimal working example is a collection of source code and other data files which allow a bug or problem to be demonstrated and reproduced. The important feature of a minimal working example is that it is as small and as simple as possible, such that it is just sufficient to demonstrate the problem, but without any additional complexity or dependencies which will make resolution harder.

Expected behavior:

The samples_contour function will throw: ValueError: A value in x_new is above the interpolation range.. Coming from the 

hs_samples = si(t_samples)

line at the end of the samples_contour function.

Additional context:

Just wanting to make sure I'm not approaching the sampling wrong here. The contours in the above example are coming from OMAE 2020 model in Vericron package. If this is the correct way to go about selecting the period_samples to feed into the sampling function, I have a solution in place that I'll make a PR out of.

Here's the solution just in case:

import scipy.interpolate as interp

def _samples_contour(t_samples, t_contour, hs_contour, bounds_error=True, fill_value=np.nan):
    """
    Get Hs points along a specified environmental contour using
    user-defined T values.

    Parameters
    ----------
    t_samples : np.array
        Points for sampling along return contour
    t_contour : np.array
        T values along contour
    hs_contour : np.array
        Hs values along contour
    bounds_error: bool
        Error on values outside of bounds
    fill_value: float = np.nan
        Fill value if bounds_error is False.  Fills in when interpolation fails for a given t_samples point.

    Returns
    -------
    hs_samples : nparray
        points sampled along return contour
    """
    assert isinstance(t_samples, np.ndarray), "t_samples must be of type np.ndarray"
    assert isinstance(t_contour, np.ndarray), "t_contour must be of type np.ndarray"
    assert isinstance(hs_contour, np.ndarray), "hs_contour must be of type np.ndarray"
    assert isinstance(bounds_error, bool), "bounds_error flag must be of type bool"
    if not bounds_error:
        try:
            fill_value = float(fill_value)
        except: 
            raise AssertionError( "if bounds_error is false, fill_value must be a float or np.nan")

    # finds minimum and maximum energy period values
    amin = np.argmin(t_contour)
    amax = np.argmax(t_contour)
    aamin = np.min([amin, amax])
    aamax = np.max([amin, amax])
    # finds points along the contour
    w1 = hs_contour[aamin:aamax]
    w2 = np.concatenate((hs_contour[aamax:], hs_contour[:aamin]))
    if np.max(w1) > np.max(w2):
        x1 = t_contour[aamin:aamax]
        y1 = hs_contour[aamin:aamax]
    else:
        x1 = np.concatenate((t_contour[aamax:], t_contour[:aamin]))
        y1 = np.concatenate((hs_contour[aamax:], hs_contour[:aamin]))
    # sorts data based on the max and min energy period values
    ms = np.argsort(x1)
    x = x1[ms]
    y = y1[ms]
    # interpolates the sorted data
    si = interp.interp1d(x, y, bounds_error=bounds_error, fill_value=fill_value)
    # finds the wave height based on the user specified energy period values
    hs_samples = si(t_samples)

    return hs_samples
ssolson commented 1 year ago

@Graham-EGI thank you for submitting this report.

I looked into it and I believe the issue is here: https://github.com/MHKiT-Software/MHKiT-Python/blob/master/mhkit/wave/contours.py#L1754

As the slicing is non-inclusive of the final element (e.g. aamax).

Can you confirm if this is the same reason you found for the error and are still planning to submit a PR?

thanks again 🙏

Graham-EGI commented 1 year ago

@ssolson thanks for the quick reply! Yea that's what I thought the issue was too.

I'm still planning to submit this PR, and will do so this week. The above solution in mind still preserves the existing behavior, just allows the user to "silence" the errors if wanting to just throw the full Period range into np.arrange() or linspace().

No need to keep the issue open, unless it's helpful for me to reference in the upcoming PR!

ssolson commented 8 months ago

@Graham-EGI I wanted to follow up on this issue. As far as I can tell this was never fixed. Do you have an unsubmitted PR or would it be better if I submitted one?

johannes-spinneken commented 8 months ago

@ssolson - Graham no longer works with EGI. Perhaps @mbruggs can comment on this?

ssolson commented 7 months ago

Resolved in #278