Plotting

modelskill.plotting

The plotting module provides functions useful for skill assessment that can be used independently of the comparison module.

• scatter is a function that can be used to plot a scatter suitable for skill assessment, with a 1:1 line and a linear regression line.
• wind_rose is a function that can be used to plot a dual wind rose to compare two datasets of magnitudes and directions.
• spatial_overview is a function that can be used to plot a spatial overview of two datasets.
• temporal_coverage is a function that can be used to plot the temporal coverage of two datasets.

scatter

scatter(x, y, *, bins=120, quantiles=None, fit_to_quantiles=False, show_points=None, show_hist=None, show_density=None, norm=None, backend='matplotlib', figsize=(8, 8), xlim=None, ylim=None, reg_method='ols', title='', xlabel='', ylabel='', skill_table=False, skill_scores=None, skill_score_unit='', ax=None, **kwargs)

Scatter plot showing compared data: observation vs modelled Optionally, with density histogram.

Parameters:

Name	Type	Description	Default
x	ndarray	X values e.g model values, must be same length as y	required
y	ndarray	Y values e.g observation values, must be same length as x	required
bins	int | float	bins for the 2D histogram on the background. By default 120 bins. if int, represents the number of bins of 2D if float, represents the bin size if sequence (list of int or float), represents the bin edges	120
quantiles	int | Sequence[float] | None	number of quantiles for QQ-plot, by default None and will depend on the scatter data length (10, 100 or 1000) if int, this is the number of points if sequence (list of floats), represents the desired quantiles (from 0 to 1)	None
fit_to_quantiles	bool	by default the regression line is fitted to all data, if True, it is fitted to the quantiles which can be useful to represent the extremes of the distribution, by default False	False
show_points	(bool, int, float)	Should the scatter points be displayed? None means: show all points if fewer than 1e4, otherwise show 1e4 sample points, by default None. float: fraction of points to show on plot from 0 to 1. eg 0.5 shows 50% of the points. int: if 'n' (int) given, then 'n' points will be displayed, randomly selected.	None
show_hist	bool	show the data density as a 2d histogram, by default None	None
show_density	Optional[bool]	show the data density as a colormap of the scatter, by default None. If both show_density and show_hist are None, then show_density is used by default. for binning the data, the previous kword bins=Float is used	None
norm	Normalize	colormap normalization If None, defaults to matplotlib.colors.PowerNorm(vmin=1,gamma=0.5)	None
backend	str	use "plotly" (interactive) or "matplotlib" backend, by default "matplotlib"	'matplotlib'
figsize	tuple	width and height of the figure, by default (8, 8)	(8, 8)
xlim	tuple	plot range for the observation (xmin, xmax), by default None	None
ylim	tuple	plot range for the model (ymin, ymax), by default None	None
reg_method	str or bool	method for determining the regression line "ols" : ordinary least squares regression "odr" : orthogonal distance regression, False : no regression line by default "ols"	'ols'
title	str	plot title, by default None	''
xlabel	str	x-label text on plot, by default None	''
ylabel	str	y-label text on plot, by default None	''
skill_table	Optional[str | Sequence[str] | bool]	calculate skill scores and show in box next to the plot, True will show default metrics, list of metrics will show these skill scores, by default False, Note: cannot be used together with skill_scores argument	False
skill_scores	dict[str, float]	dictionary with skill scores to be shown in box next to the plot, by default None Note: cannot be used together with skill_table argument	None
skill_score_unit	str	unit for skill_scores, by default None	''
ax	Axes	axes to plot on, by default None	None
**kwargs			{}

Returns:

Type	Description
Axes	The axes on which the scatter plot was drawn.

Source code in modelskill/plotting/_scatter.py

def scatter(
    x: np.ndarray,
    y: np.ndarray,
    *,
    bins: int | float = 120,
    quantiles: int | Sequence[float] | None = None,
    fit_to_quantiles: bool = False,
    show_points: bool | int | float | None = None,
    show_hist: Optional[bool] = None,
    show_density: Optional[bool] = None,
    norm: Optional[colors.Normalize] = None,
    backend: Literal["matplotlib", "plotly"] = "matplotlib",
    figsize: Tuple[float, float] = (8, 8),
    xlim: Optional[Tuple[float, float]] = None,
    ylim: Optional[Tuple[float, float]] = None,
    reg_method: str | bool = "ols",
    title: str = "",
    xlabel: str = "",
    ylabel: str = "",
    skill_table: Optional[str | Sequence[str] | bool] = False,
    skill_scores: Mapping[str, float] | None = None,
    skill_score_unit: Optional[str] = "",
    ax: Optional[Axes] = None,
    **kwargs,
) -> Axes:
    """Scatter plot showing compared data: observation vs modelled
    Optionally, with density histogram.

    Parameters
    ----------
    x: np.array
        X values e.g model values, must be same length as y
    y: np.array
        Y values e.g observation values, must be same length as x
    bins: (int, float, sequence), optional
        bins for the 2D histogram on the background. By default 120 bins.
        if int, represents the number of bins of 2D
        if float, represents the bin size
        if sequence (list of int or float), represents the bin edges
    quantiles: (int, sequence), optional
        number of quantiles for QQ-plot, by default None and will depend on the scatter data length (10, 100 or 1000)
        if int, this is the number of points
        if sequence (list of floats), represents the desired quantiles (from 0 to 1)
    fit_to_quantiles: bool, optional
        by default the regression line is fitted to all data, if True, it is fitted to the quantiles
        which can be useful to represent the extremes of the distribution, by default False
    show_points : (bool, int, float), optional
        Should the scatter points be displayed?
        None means: show all points if fewer than 1e4, otherwise show 1e4 sample points, by default None.
        float: fraction of points to show on plot from 0 to 1. eg 0.5 shows 50% of the points.
        int: if 'n' (int) given, then 'n' points will be displayed, randomly selected.
    show_hist : bool, optional
        show the data density as a 2d histogram, by default None
    show_density: bool, optional
        show the data density as a colormap of the scatter, by default None. If both `show_density` and `show_hist`
        are None, then `show_density` is used by default.
        for binning the data, the previous kword `bins=Float` is used
    norm : matplotlib.colors.Normalize
        colormap normalization
        If None, defaults to matplotlib.colors.PowerNorm(vmin=1,gamma=0.5)
    backend : str, optional
        use "plotly" (interactive) or "matplotlib" backend, by default "matplotlib"
    figsize : tuple, optional
        width and height of the figure, by default (8, 8)
    xlim : tuple, optional
        plot range for the observation (xmin, xmax), by default None
    ylim : tuple, optional
        plot range for the model (ymin, ymax), by default None
    reg_method : str or bool, optional
        method for determining the regression line
        "ols" : ordinary least squares regression
        "odr" : orthogonal distance regression,
        False : no regression line
        by default "ols"
    title : str, optional
        plot title, by default None
    xlabel : str, optional
        x-label text on plot, by default None
    ylabel : str, optional
        y-label text on plot, by default None
    skill_table: str, List[str], bool, optional
        calculate skill scores and show in box next to the plot,
        True will show default metrics, list of metrics will show
        these skill scores, by default False,
        Note: cannot be used together with skill_scores argument
    skill_scores : dict[str, float], optional
        dictionary with skill scores to be shown in box next to
        the plot, by default None
        Note: cannot be used together with skill_table argument
    skill_score_unit : str, optional
        unit for skill_scores, by default None
    ax : matplotlib.axes.Axes, optional
        axes to plot on, by default None
    **kwargs

    Returns
    -------
    matplotlib.axes.Axes
        The axes on which the scatter plot was drawn.
    """
    if "skill_df" in kwargs:
        warnings.warn(
            "The `skill_df` keyword argument is deprecated. Use `skill_scores` instead.",
            FutureWarning,
        )
        skill_scores = kwargs.pop("skill_df").to_dict("records")[0]

    if show_hist is None and show_density is None:
        # Default: points density
        show_density = True

    if len(x) != len(y):
        raise ValueError("x & y are not of equal length")

    if norm is None:
        norm = colors.PowerNorm(vmin=1, gamma=0.5)

    x_sample, y_sample = sample_points(x, y, show_points)
    xq, yq = quantiles_xy(x, y, quantiles)

    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()
    xymin = min([xmin, ymin])
    xymax = max([xmax, ymax])

    nbins_hist, binsize = _get_bins(bins, xymin=xymin, xymax=xymax)

    if xlim is None:
        xlim = (xymin - binsize, xymax + binsize)

    if ylim is None:
        ylim = (xymin - binsize, xymax + binsize)

    x_trend = np.array([xlim[0], xlim[1]])

    if show_hist and show_density:
        raise TypeError(
            "if `show_hist=True` then `show_density` must be either `False` or `None`"
        )

    z = None
    if show_density and len(x_sample) > 0:
        if not isinstance(bins, (float, int)):
            raise TypeError(
                "if `show_density=True` then bins must be either float or int"
            )

        # calculate density data
        z = __scatter_density(x_sample, y_sample, binsize=binsize)
        idx = z.argsort()
        # Sort data by colormaps
        x_sample, y_sample, z = x_sample[idx], y_sample[idx], z[idx]
        # scale Z by sample size
        z = z * len(x) / len(x_sample)

    PLOTTING_BACKENDS: dict[str, Callable] = {
        "matplotlib": _scatter_matplotlib,
        "plotly": _scatter_plotly,
    }

    if backend not in PLOTTING_BACKENDS:
        raise ValueError(f"backend must be one of {list(PLOTTING_BACKENDS.keys())}")

    if skill_table:
        from modelskill import from_matched

        if skill_scores is not None:
            raise ValueError(
                "Cannot pass skill_scores and skill_table at the same time"
            )
        df = pd.DataFrame({"obs": x, "model": y})
        cmp = from_matched(df)
        metrics = None if skill_table is True else skill_table
        skill = cmp.skill(metrics=metrics)
        skill_scores = skill.to_dict("records")[0]

    return PLOTTING_BACKENDS[backend](
        x=x,
        y=y,
        x_sample=x_sample,
        y_sample=y_sample,
        z=z,
        xq=xq,
        yq=yq,
        x_trend=x_trend,
        show_density=show_density,
        norm=norm,
        show_points=show_points,
        show_hist=show_hist,
        nbins_hist=nbins_hist,
        reg_method=reg_method,
        xlabel=xlabel,
        ylabel=ylabel,
        figsize=figsize,
        xlim=xlim,
        ylim=ylim,
        title=title,
        skill_scores=skill_scores,
        skill_score_unit=skill_score_unit,
        fit_to_quantiles=fit_to_quantiles,
        ax=ax,
        **kwargs,
    )

spatial_overview

spatial_overview(obs, mod=None, ax=None, figsize=None, title=None)

Plot observation points on a map showing the model domain

Parameters:

Name	Type	Description	Default
obs	List[Observation]	List of observations to be shown on map	required
mod	Union[ModelResult, GeometryFM]	Model domain to be shown as outline	None
ax		Adding to existing axis, instead of creating new fig	None
figsize	(float, float)	figure size, by default None	None
title	Optional[str]	plot title, default empty	None

See Also

temporal_coverage

Returns:

Type	Description
Axes	The matplotlib axes object

Examples:

>>> import modelskill as ms
>>> o1 = ms.PointObservation('HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
>>> o2 = ms.TrackObservation("Alti_c2_Dutch.dfs0", item=3, name="c2")
>>> mr1 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
>>> mr2 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)
>>> ms.plotting.spatial_overview([o1, o2], [mr1, mr2])

Source code in modelskill/plotting/_spatial_overview.py

def spatial_overview(
    obs: List[Observation],
    mod=None,
    ax=None,
    figsize: Optional[Tuple] = None,
    title: Optional[str] = None,
) -> matplotlib.axes.Axes:
    """Plot observation points on a map showing the model domain

    Parameters
    ----------
    obs: list[Observation]
        List of observations to be shown on map
    mod : Union[ModelResult, mikeio.GeometryFM], optional
        Model domain to be shown as outline
    ax: matplotlib.axes, optional
        Adding to existing axis, instead of creating new fig
    figsize : (float, float), optional
        figure size, by default None
    title: str, optional
        plot title, default empty

    See Also
    --------
    temporal_coverage

    Returns
    -------
    matplotlib.axes.Axes
        The matplotlib axes object

    Examples
    --------
    >>> import modelskill as ms
    >>> o1 = ms.PointObservation('HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
    >>> o2 = ms.TrackObservation("Alti_c2_Dutch.dfs0", item=3, name="c2")
    >>> mr1 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
    >>> mr2 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)
    >>> ms.plotting.spatial_overview([o1, o2], [mr1, mr2])
    """
    obs = [] if obs is None else list(obs) if isinstance(obs, Sequence) else [obs]  # type: ignore
    mod = [] if mod is None else list(mod) if isinstance(mod, Sequence) else [mod]  # type: ignore

    ax = _get_ax(ax=ax, figsize=figsize)
    offset_x = 1  # TODO: better default

    for m in mod:
        # TODO: support Gridded ModelResults
        if isinstance(m, (PointModelResult, TrackModelResult)):
            raise ValueError(
                f"Model type {type(m)} not supported. Only DfsuModelResult and mikeio.GeometryFM supported!"
            )
        if hasattr(m, "data") and hasattr(m.data, "geometry"):
            # mod_name = m.name  # TODO: better support for multiple models
            m = m.data.geometry
        if hasattr(m, "node_coordinates"):
            xn = m.node_coordinates[:, 0]
            offset_x = 0.02 * (max(xn) - min(xn))
        m.plot.outline(ax=ax)

    for o in obs:
        if isinstance(o, PointObservation):
            ax.scatter(x=o.x, y=o.y, marker="x")
            ax.annotate(o.name, (o.x + offset_x, o.y))  # type: ignore
        elif isinstance(o, TrackObservation):
            if o.n_points < 10000:
                ax.scatter(x=o.x, y=o.y, c=o.values, marker=".", cmap="Reds")
            else:
                print(f"{o.name}: Too many points to plot")
                # TODO: group by lonlat bin or sample randomly
        else:
            raise ValueError(
                f"Could not show observation {o}. Only PointObservation and TrackObservation supported."
            )

    if not title:
        title = "Spatial coverage"
    ax.set_title(title)

    return ax

taylor_diagram

taylor_diagram(obs_std, points, figsize=(7, 7), obs_text='Observations', normalize_std=False, ax=None, title='Taylor diagram')

Plot a Taylor diagram using the given observations and points.

Parameters:

Name	Type	Description	Default
obs_std	float	Standard deviation of the observations.	required
points	list of TaylorPoint objects or a single TaylorPoint object	Points to plot on the Taylor diagram.	required
figsize	tuple	Figure size in inches. Default is (7, 7).	(7, 7)
obs_text	str	Label for the observations. Default is "Observations".	'Observations'
normalize_std	bool	Whether to normalize the standard deviation of the points by the standard deviation of the observations. Default is False.	False
title	str	Title of the plot. Default is "Taylor diagram".	'Taylor diagram'

Returns:

Type	Description
Figure	The matplotlib figure object

Source code in modelskill/plotting/_taylor_diagram.py

def taylor_diagram(
    obs_std,
    points,
    figsize=(7, 7),
    obs_text="Observations",
    normalize_std=False,
    ax=None,
    title="Taylor diagram",
) -> matplotlib.figure.Figure:
    """
    Plot a Taylor diagram using the given observations and points.

    Parameters
    -----------
    obs_std : float
        Standard deviation of the observations.
    points : list of TaylorPoint objects or a single TaylorPoint object
        Points to plot on the Taylor diagram.
    figsize : tuple, optional
        Figure size in inches. Default is (7, 7).
    obs_text : str, optional
        Label for the observations. Default is "Observations".
    normalize_std : bool, optional
        Whether to normalize the standard deviation of the points by the standard deviation of the observations. Default is False.
    title : str, optional
        Title of the plot. Default is "Taylor diagram".

    Returns
    --------
    matplotlib.figure.Figure
            The matplotlib figure object
    """

    if np.isscalar(figsize):
        figsize = (figsize, figsize)
    elif figsize[0] != figsize[1]:
        warnings.warn(
            "It is strongly recommended that the aspect ratio is 1:1 for Taylor diagrams"
        )
    fig = plt.figure(figsize=figsize)

    # srange=(0, 1.5),
    if len(obs_text) > 30:
        obs_text = obs_text[:25] + "..."

    td = TaylorDiagram(
        obs_std, fig=fig, rect=111, label=obs_text, normalize_std=normalize_std
    )
    contours = td.add_contours(levels=8, colors="0.5", linestyles="dotted")
    plt.clabel(contours, inline=1, fontsize=10, fmt="%.2f")

    if isinstance(points, TaylorPoint):
        points = [points]
    for p in points:
        assert isinstance(p, TaylorPoint)
        m = "o" if p.marker is None else p.marker
        ms = "6" if p.marker_size is None else p.marker_size
        std = p.std / p.obs_std if normalize_std else p.std
        td.add_sample(std, p.cc, marker=m, ms=ms, ls="", label=p.name)
        # marker=f"${1}$",
        # td.add_sample(0.2, 0.8, marker="+", ms=15, mew=1.2, ls="", label="m2")
    td.add_grid()
    fig.legend(
        td.samplePoints,
        [p.get_label() for p in td.samplePoints],
        numpoints=1,
        prop=dict(size="medium"),
        loc="upper right",
    )
    fig.suptitle(title, size="x-large")

    # prevent the plot from being displayed, since it is also displayed by the returned object
    plt.close()
    return fig

temporal_coverage

temporal_coverage(obs=None, mod=None, *, limit_to_model_period=True, marker='_', ax=None, figsize=None, title=None)

Plot graph showing temporal coverage for all observations and models

Parameters:

Name	Type	Description	Default
obs	List[Observation]	Show observation(s) as separate lines on plot	None
mod	List[ModelResult]	Show model(s) as separate lines on plot, by default None	None
limit_to_model_period	bool	Show temporal coverage only for period covered by the model, by default True	True
marker	str	plot marker for observations, by default "_"	'_'
ax		Adding to existing axis, instead of creating new fig	None
figsize	Tuple(float, float)	size of figure, by default (7, 0.45*n_lines)	None
title		plot title, default empty	None

See Also

spatial_overview

Returns:

Type	Description
Axes	The matplotlib axes object

Examples:

>>> import modelskill as ms
>>> o1 = ms.PointObservation('HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
>>> o2 = ms.TrackObservation("Alti_c2_Dutch.dfs0", item=3, name="c2")
>>> mr1 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
>>> mr2 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)
>>> ms.plotting.temporal_coverage([o1, o2], [mr1, mr2])
>>> ms.plotting.temporal_coverage([o1, o2], mr2, limit_to_model_period=False)
>>> ms.plotting.temporal_coverage(o2, [mr1, mr2], marker=".")
>>> ms.plotting.temporal_coverage(mod=[mr1, mr2], figsize=(5,3))

Source code in modelskill/plotting/_temporal_coverage.py

def temporal_coverage(
    obs=None,
    mod=None,
    *,
    limit_to_model_period=True,
    marker="_",
    ax=None,
    figsize=None,
    title=None,
) -> matplotlib.axes.Axes:
    """Plot graph showing temporal coverage for all observations and models

    Parameters
    ----------
    obs : List[Observation], optional
        Show observation(s) as separate lines on plot
    mod : List[ModelResult], optional
        Show model(s) as separate lines on plot, by default None
    limit_to_model_period : bool, optional
        Show temporal coverage only for period covered
        by the model, by default True
    marker : str, optional
        plot marker for observations, by default "_"
    ax: matplotlib.axes, optional
        Adding to existing axis, instead of creating new fig
    figsize : Tuple(float, float), optional
        size of figure, by default (7, 0.45*n_lines)
    title: str, optional
        plot title, default empty

    See Also
    --------
    spatial_overview

    Returns
    -------
    matplotlib.axes.Axes
        The matplotlib axes object

    Examples
    --------
    >>> import modelskill as ms
    >>> o1 = ms.PointObservation('HKNA_Hm0.dfs0', item=0, x=4.2420, y=52.6887, name="HKNA")
    >>> o2 = ms.TrackObservation("Alti_c2_Dutch.dfs0", item=3, name="c2")
    >>> mr1 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast.dfsu', name='SW_1', item=0)
    >>> mr2 = ms.DfsuModelResult('HKZN_local_2017_DutchCoast_v2.dfsu', name='SW_2', item=0)
    >>> ms.plotting.temporal_coverage([o1, o2], [mr1, mr2])
    >>> ms.plotting.temporal_coverage([o1, o2], mr2, limit_to_model_period=False)
    >>> ms.plotting.temporal_coverage(o2, [mr1, mr2], marker=".")
    >>> ms.plotting.temporal_coverage(mod=[mr1, mr2], figsize=(5,3))
    """
    obs = [] if obs is None else list(obs) if isinstance(obs, Sequence) else [obs]
    mod = [] if mod is None else list(mod) if isinstance(mod, Sequence) else [mod]

    n_lines = len(obs) + len(mod)
    if figsize is None:
        ysize = max(2.0, 0.45 * n_lines)
        figsize = (7, ysize)

    fig, ax = _get_fig_ax(ax=ax, figsize=figsize)
    y = np.repeat(0.0, 2)
    labels = []

    if len(mod) > 0:
        for mr in mod:
            y += 1.0
            plt.plot([mr.time[0], mr.time[-1]], y)
            labels.append(mr.name)

    for o in obs:
        y += 1.0
        plt.plot(o.time, y[0] * np.ones(len(o.time)), marker, markersize=5)
        labels.append(o.name)

    if len(mod) > 0 and limit_to_model_period:
        mr = mod[0]  # take first model
        plt.xlim([mr.time[0], mr.time[-1]])

    plt.yticks(np.arange(n_lines) + 1, labels)
    if len(mod) > 0:
        for j in range(len(mod)):
            ax.get_yticklabels()[j].set_fontstyle("italic")
            ax.get_yticklabels()[j].set_weight("bold")
            # set_color("#004165")
    fig.autofmt_xdate()

    if title:
        ax.set_title(title)
    return ax

wind_rose

wind_rose(data, *, labels=('Measurement', 'Model'), mag_step=None, n_sectors=16, calm_threshold=None, calm_size=None, calm_text='Calm', r_step=0.1, r_max=None, legend=True, cmap1='viridis', cmap2='Greys', mag_bins=None, max_bin=None, n_dir_labels=None, secondary_dir_step_factor=2.0, figsize=(8, 8), ax=None, title=None)

Plots a (dual) wind (wave or current) roses with calms.

The size of the calm is determined by the primary (measurement) data.

Parameters:

Name	Type	Description	Default
data		array with 2 or 4 columns (magnitude, direction, magnitude2, direction2)	required
labels		labels for the legend(s)	('Measurement', 'Model')
mag_step	Optional[float]	discretization for magnitude (delta_r, in radial direction )	None
n_sectors	int	number of directional sectors	16
calm_threshold	Optional[float]	minimum value for data being counted as valid (i.e. below this is calm)	None
calm_text	str	text to display in calm.	'Calm'
r_step	float	radial axis discretization. By default 0.1 i.e. every 10%.	0.1
r_max	Optional[float]	maximum radius (%) of plot, e.g. if 50% wanted then r_max=0.5	None
max_bin	Optional[float]	max value to truncate the data, e.g., max_bin=1.0 if hm0=1m is the desired final bin.	None
mag_bins	array of floats (optional) Default = None	force bins to array of values, e.g. when specifying non-equidistant bins.	None
legend	bool	show legend	True
cmap1	string. Default= 'viridis'	colormap for main axis	'viridis'
cmap2	string. Default= 'Greys'	colormap for secondary axis	'Greys'
n_dir_labels	int. Default= 4	number of labels in the polar plot, choose between 4, 8 or 16, default is to use the same as n_sectors	None
secondary_dir_step_factor	float. Default= 2.0	reduce width of secondary axis by this factor	2.0
figsize	Tuple[float, float]	figure size	(8, 8)
ax		Matplotlib axis to plot on defined as polar, it can be done using "subplot_kw = dict(projection = 'polar')". Default = None, new axis created.	None
title		title of the plot	None

Returns:

Type	Description
Axes	Matplotlib axis with the plot

Source code in modelskill/plotting/_wind_rose.py

def wind_rose(
    data,
    *,
    labels=("Measurement", "Model"),
    mag_step: Optional[float] = None,
    n_sectors: int = 16,
    calm_threshold: Optional[float] = None,  # TODO rename to vmin?
    calm_size: Optional[float] = None,
    calm_text: str = "Calm",
    r_step: float = 0.1,
    r_max: Optional[float] = None,
    legend: bool = True,
    cmap1: str = "viridis",
    cmap2: str = "Greys",
    mag_bins: Optional[List[float]] = None,
    max_bin: Optional[float] = None,  # TODO rename to vmax?
    n_dir_labels: Optional[int] = None,
    secondary_dir_step_factor: float = 2.0,
    figsize: Tuple[float, float] = (8, 8),
    ax=None,
    title=None,
) -> matplotlib.axes.Axes:
    """Plots a (dual) wind (wave or current) roses with calms.

    The size of the calm is determined by the primary (measurement) data.

    Parameters
    ----------
    data: array-like
        array with 2 or 4 columns (magnitude, direction, magnitude2, direction2)
    labels: tuple of strings. Default= ("Measurement", "Model")
        labels for the legend(s)
    mag_step: float, (optional) Default= None
        discretization for magnitude (delta_r, in radial direction )
    n_sectors: int (optional) Default= 16
        number of directional sectors
    calm_threshold: float (optional) Default= None (auto calculated)
        minimum value for data being counted as valid (i.e. below this is calm)
    calm_text: str (optional) Default: 'Calm'
        text to display in calm.
    r_step: float (optional) Default= 0.1
        radial axis discretization. By default 0.1 i.e. every 10%.
    r_max: float (optional) Default= None
        maximum radius (%) of plot, e.g. if 50% wanted then r_max=0.5
    max_bin:  float (optional) Default= None
        max value to truncate the data, e.g.,  max_bin=1.0 if hm0=1m is the desired final bin.
    mag_bins : array of floats (optional) Default = None
        force bins to array of values, e.g. when specifying non-equidistant bins.
    legend: boolean. Default= True
        show legend
    cmap1 : string. Default= 'viridis'
        colormap for main axis
    cmap2 : string. Default= 'Greys'
        colormap for secondary axis
    n_dir_labels : int. Default= 4
        number of labels in the polar plot, choose between 4, 8 or 16, default is to use the same as n_sectors
    secondary_dir_step_factor : float. Default= 2.0
        reduce width of secondary axis by this factor
    figsize: tuple(float,float)
        figure size
    ax: Matplotlib axis Default= None
        Matplotlib axis to plot on defined as polar, it can be done using "subplot_kw = dict(projection = 'polar')". Default = None, new axis created.
    title: str Default= None
        title of the plot

    Returns
    -------
    matplotlib.axes.Axes
        Matplotlib axis with the plot
    """
    if hasattr(data, "to_numpy"):
        data = data.to_numpy()

    # check that data is array_like
    assert hasattr(data, "__array__"), "data must be array_like"

    data_1 = data[:, 0:2]  # primary magnitude and direction
    magmax = data_1[:, 0].max()

    ncols = data.shape[1]
    assert ncols in [2, 4], "data must have 2 or 4 columns"
    dual = ncols == 4

    if dual:
        data_2 = data[:, 2:4]  # secondary magnitude and direction
        magmax = max(magmax, data_2[:, 0].max())
        assert len(labels) == 2, "labels must have 2 elements"

    # magnitude bins
    ui, vmin, vmax = pretty_intervals(
        magmax,
        mag_bins,
        mag_step,
        calm_threshold,
        max_bin,
    )

    dir_step = 360 // n_sectors

    if n_dir_labels is None:
        if n_sectors in (4, 8, 16):
            n_dir_labels = n_sectors
        else:
            # Directional labels are not identical to the number of sectors, use a sane default
            n_dir_labels = 16

    dh = _dirhist2d(data_1, ui=ui, dir_step=dir_step)
    calm = dh.calm

    if dual:
        assert len(data_1) == len(data_2), "data_1 and data_2 must have same length"
        dh2 = _dirhist2d(data_2, ui=ui, dir_step=dir_step)
        assert dh.density.shape == dh2.density.shape

    ri, rmax = _calc_radial_ticks(counts=dh.density, step=r_step, stop=r_max)

    # Resize calm
    # TODO this overwrites the calm value calculated above
    if calm_size is not None:
        calm = calm_size

    cmap = _get_cmap(cmap1)

    if ax is None:
        _, ax = plt.subplots(figsize=figsize, subplot_kw=dict(projection="polar"))

    ax.set_title(title)
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)

    dir_labels = directional_labels(n_dir_labels)
    grid = np.linspace(0, 360, n_dir_labels + 1)[:-1]
    ax.set_thetagrids(grid, dir_labels)

    # ax.tick_params(pad=-24)

    ax.set_ylim(0, calm + rmax)
    ax.set_yticks(ri + calm)
    tick_labels = [f"{tick * 100 :.0f}%" for tick in ri]
    ax.set_yticklabels(tick_labels)
    ax.set_rlabel_position(5)

    if vmin > 0:
        _add_calms_to_ax(ax, threshold=calm, text=calm_text)

    # primary histogram (model)
    p = _create_patch(
        thetac=dh.dir_centers,
        dir_step=dir_step,
        calm=calm,
        ui=ui,
        counts=dh.density,
        cmap=cmap,
        vmax=vmax,
    )
    ax.add_collection(p)

    if legend:
        _add_legend_to_ax(
            ax,
            cmap=cmap,
            vmax=vmax,
            ui=ui,
            calm=calm,
            counts=dh.density,
            label=labels[0],
            primary=True,
            dual=dual,
        )

    if dual:
        # add second histogram (observation)
        cmap = _get_cmap(cmap2)

        # TODO should this be calm2?
        p = _create_patch(
            thetac=dh.dir_centers,
            dir_step=dir_step,
            calm=calm,
            ui=ui,
            counts=dh2.density,
            cmap=cmap,
            vmax=vmax,
            dir_step_factor=secondary_dir_step_factor,
        )
        ax.add_collection(p)

        if legend:
            _add_legend_to_ax(
                ax,
                cmap=cmap,
                vmax=vmax,
                ui=ui,
                calm=dh2.calm,
                counts=dh2.density,
                label=labels[1],
                primary=False,
                dual=dual,
            )

    return ax