Wind speed estimate at hub-height using shear

[elenya-grant]

Wind speed estimate at hub-height using shear

Add function for estimating wind resource at the turbine hub-height using shear. This should be included as a new option that is specific to wind speed (but similar to resource_data_averaging_method) in the FlorisWindPlantPerformanceConfig.

Proposed solution

Theres a few possible options/use-cases:

1. wind speed only available at one resource height, user must provide wind shear value (psi)

ws_est = uref*(hub_height/zref)**psi

2. wind speed available at multiple heights, use the resource heights nearest to the hub-height to estimate shear, then using this shear value, estimate the wind speed at the hub-height

def estimate_wind_speed(u0, z0, uref, zref, hub_height):
    psi = np.log(z0/zref)/np.log(u0/uref)
    ws_est = uref*(hub_height/zref)**psi
    return ws_est

3. wind speed available at multiple heights, use the resource heights nearest to the hub-height to estimate curve coefficients for wind speed from hub-height and use this to estimate wind speed at the hub-height. This can be done in two different ways:
   • one curve-fit for the entire year
   • one curve fit per time-step, although this may require using more than just the bounding heights.

def height_to_winspeed_func(ur_zr_z, psi):
    ur, zr, z = ur_zr_z
    u = ur * (z / zr) ** psi
    return u

For the entire year:

def simple_estimate_wind_speed_with_curve_fit(
    wind_resource_data: dict,
    bounding_resource_heights: tuple[int] | list[int],
    hub_height: float | int,
    ):
    ws_heights = np.array(bounding_resource_heights)
    n_timesteps = len(wind_resource_data[f"wind_speed_{int(ws_heights[0])}m"])
    # calc closest height
    ub_diff = np.abs(np.max(ws_heights) - hub_height)
    lb_diff = np.abs(np.min(ws_heights) - hub_height)

    if ub_diff >= lb_diff:
        # lower-bound is closer, use lower bound as reference and upper bound as input
        z_ref = np.min(ws_heights) * np.ones(n_timesteps)
        ws_ref = wind_resource_data[f"wind_speed_{int(np.min(ws_heights))}m"]
        z = np.max(ws_heights) * np.ones(n_timesteps)
        ws = wind_resource_data[f"wind_speed_{int(np.max(ws_heights))}m"]

    else:
        # upper bound is closer, use upper bound as reference and lower bound as input
        z_ref = np.max(ws_heights) * np.ones(n_timesteps)
        ws_ref = wind_resource_data[f"wind_speed_{int(np.max(ws_heights))}m"]
        z = np.min(ws_heights) * np.ones(n_timesteps)
        ws = wind_resource_data[f"wind_speed_{int(np.min(ws_heights))}m"]

    curve_coeff, curve_cov = scipy.optimize.curve_fit(
        height_to_winspeed_func,
        (ws_ref, z_ref, z),
        ws,
        p0=(1.0),
    )
    ws_at_hubheight = height_to_winspeed_func(
        (ws_ref, z_ref, hub_height * np.ones(n_timesteps)), *curve_coeff
    )
    return ws_at_hubheight

Estimated per timestep, although this code does not work at the moment:

def estimate_wind_speed_with_curve_fit(
    wind_resource_data: dict,
    bounding_resource_heights: tuple[int] | list[int],
    hub_height: float | int,
):
    ws_dict = {k: v for k, v in wind_resource_data.items() if "wind_speed" in k}
    ws_heights = np.array(bounding_resource_heights)
    n_timesteps = len(ws_dict[f"wind_speed_{int(ws_heights[0])}m"])

    # calc closest height
    ub_diff = np.abs(np.max(ws_heights) - hub_height)
    lb_diff = np.abs(np.min(ws_heights) - hub_height)

    if ub_diff >= lb_diff:
        # lower-bound is closer, use lower bound as reference and upper bound as input
        z_ref = np.min(ws_heights) * np.ones(n_timesteps)
        ws_ref = ws_dict[f"wind_speed_{int(np.min(ws_heights))}m"]
        z = np.max(ws_heights) * np.ones(n_timesteps)
        ws = ws_dict[f"wind_speed_{int(np.max(ws_heights))}m"]

    else:
        # upper bound is closer, use upper bound as reference and lower bound as input
        z_ref = np.max(ws_heights) * np.ones(n_timesteps)
        ws_ref = ws_dict[f"wind_speed_{int(np.max(ws_heights))}m"]
        z = np.min(ws_heights) * np.ones(n_timesteps)
        ws = ws_dict[f"wind_speed_{int(np.min(ws_heights))}m"]
    
ws_at_hubheight = np.zeros(n_timesteps)
    for i in range(n_timesteps):
        # error is thrown here for some reason
        curve_coeff, curve_cov = scipy.optimize.curve_fit(
            height_to_winspeed_func,
            (np.array(ws_ref[i]), np.array(z_ref[i]), np.array(z[i])),
            np.array(ws[i]),
            p0=(1.0),
        )
        ws_at_hubheight[i] = height_to_winspeed_func(
            (ws_ref[i], z_ref[i], hub_height), *curve_coeff
        )

Alternatives considered

Additional context

Brought up with the introduction of h2integrate/converters/wind/tools/resource_tools.py in PR #372