Speed up Solar tilt_orientation_factor plot method

[Tom-Kingstone]

Description:

The fast version originally sat in #161's branch, but it should be viewed separately to that issue

def tilt_orientation_factor(
    epw_file: Path,
    ax: plt.Axes = None,
    rad_type: str = "total",
    analysis_period: AnalysisPeriod = AnalysisPeriod(),
    cmap: str = "YlOrRd",
    directions: int = 36,
    tilts: int = 9,
    quantiles: tuple[float] = (0.05, 0.25, 0.5, 0.75, 0.95),
    lims: tuple[float, float] = None,
) -> plt.Axes:
    """Create a tilt-orientation factor plot.

    Args:
        epw_file (Path): 
            The EPW file representing the weather data/location to be visualised.
        ax (plt.Axes, optional): 
            The axes to plot on. 
            Defaults to None.
        rad_type (str, optional): 
            The type of radiation to plot. 
            Defaults to "total", with options of "total", "direct" and "diffuse".
        analysis_period (AnalysisPeriod, optional): 
            The analysis period over which radiation shall be summarised. 
            Defaults to AnalysisPeriod().
        cmap (str, optional): 
            The colormap to apply. 
            Defaults to "YlOrRd".
        directions (int, optional): 
            The number of directions to bin data into. 
            Defaults to 36.
        tilts (int, optional): 
            The number of tilts to calculate. 
            Defaults to 9.
        quantiles (tuple[float], optional):
            The quantiles to plot. 
            Defaults to (0.05, 0.25, 0.5, 0.75, 0.95).
        lims (tuple[float, float], optional):
            The limits of the plot. 
            Defaults to None.

    Returns:
        plt.Axes: 
            The matplotlib axes.
    """

    # create dir for cached results
    _dir = Path(hb_folders.default_simulation_folder) / "_lbt_tk_solar"
    _dir.mkdir(exist_ok=True, parents=True)
    ndir = directions

    if ax is None:
        ax = plt.gca()

    cmap = plt.get_cmap(cmap)

    # create sky matrix
    smx = SkyMatrix.from_epw(
        epw_file=epw_file, high_density=True, hoys=analysis_period.hoys
    )

    # create roses per tilt angle
    _directions = np.linspace(0, 360, directions + 1)[:-1].tolist()
    _tilts = np.linspace(0, 90, tilts)[:-1].tolist() + [89.999]
    rrs: list[RadiationRose] = []
    for ta in tqdm(_tilts):
        sp = _dir / f"{epw_file.stem}_{ndir}_{ta:0.4f}.pickle"
        if sp.exists():
            rr = pickle.load(open(sp, "rb"))
        else:
            rr = RadiationRose(sky_matrix=smx, direction_count=directions, tilt_angle=ta)
            pickle.dump(rr, open(sp, "wb"))
        rrs.append(rr)
    _directions.append(360)

    # create matrix of values from results
    values = np.array([getattr(i, f"{rad_type}_values") for i in rrs])

    # repeat first result at end to close the circle
    values = values.T.tolist()
    values.append(values[0])
    values = np.array(values).T

    # create x, y coordinates per result value
    __directions, __tilts = np.meshgrid(_directions, _tilts)

    # get location of max
    _max = values.flatten().max()
    if _max == 0:
        raise ValueError(f"No solar radiation within {analysis_period}.")

    _max_idx = values.flatten().argmax()
    _max_alt = __tilts.flatten()[_max_idx]
    _max_az = __directions.flatten()[_max_idx]

    # create colormap
    if lims is None:
        norm = Normalize(vmin=0, vmax=_max)
    else:
        norm = Normalize(vmin=lims[0], vmax=lims[1])

    # create triangulation
    tri = Triangulation(x=__directions.flatten(), y=__tilts.flatten())

    # create quantile lines
    quantiles = [0.05, 0.25, 0.5, 0.75, 0.95]
    levels = [np.quantile(a=values.flatten(), q=i) for i in quantiles]
    quant_colors = [cmap(i) for i in [norm(v) for v in levels]]
    quant_colors_inv = [contrasting_color(i) for i in quant_colors]
    max_color_inv = contrasting_color(cmap(norm(_max)))

    # plot data
    tcf = ax.tricontourf(tri, values.flatten(), levels=100, cmap=cmap, norm=norm)
    tcl = ax.tricontour(
        tri,
        values.flatten(),
        levels=levels,
        colors=quant_colors_inv,
        linestyles=":",
        alpha=0.5,
    )

    # add contour labels
    def cl_fmt(x):
        return f"{x:,.0f}W/m$^2$"

    _ = ax.clabel(tcl, fontsize="small", fmt=cl_fmt)

    # add colorbar
    cb = plt.colorbar(
        tcf,
        ax=ax,
        orientation="vertical",
        drawedges=False,
        fraction=0.05,
        aspect=25,
        pad=0.02,
        label="Cumulative irradiance (W/m$^2$)",
    )
    cb.outline.set_visible(False)
    for quantile_val in levels:
        cb.ax.plot([0, 1], [quantile_val] * 2, color="k", ls="-", alpha=0.5)

    # add max-location
    ax.scatter(_max_az, _max_alt, c=max_color_inv, s=10, marker="x")
    alt_offset = (90 / 100) * 0.5 if _max_alt <= 45 else -(90 / 100) * 0.5
    az_offset = (360 / 100) * 0.5 if _max_az <= 180 else -(360 / 100) * 0.5
    ha = "left" if _max_az <= 180 else "right"
    va = "bottom" if _max_alt <= 45 else "top"
    ax.text(
        _max_az + az_offset,
        _max_alt + alt_offset,
        f"{_max:,.0f}W/m$^2$\n({_max_az:0.0f}°, {_max_alt:0.0f}°)",
        ha=ha,
        va=va,
        c=max_color_inv,
        weight="bold",
        size="small",
    )

    ax.set_xlim(0, 360)
    ax.set_ylim(0, 90)
    ax.xaxis.set_major_locator(MultipleLocator(base=30))
    ax.yaxis.set_major_locator(MultipleLocator(base=10))
    ax.set_xlabel("Panel orientation (clockwise from North at 0°)")
    ax.set_ylabel("Panel tilt (0° facing the horizon, 90° facing the sky)")

    ax.set_title(
        f"{epw_file.name}\n{rad_type.title()} irradiance (cumulative)\n{describe_analysis_period(analysis_period)}"
    )

    return ax

[Tom-Kingstone]

In addition, could apply principles from #217 and return useful information from the dataset - Could use generic collection_metadata (as the parameters for the solar tilt plot are the irradiance collections)