def threshold_brier_score(observations, forecasts, threshold, issorted=False,
                          axis=-1):
    observations = np.asarray(observations)
    threshold = np.asarray(threshold)
    forecasts = np.asarray(forecasts)
    if axis != -1:
        forecasts = move_axis_to_end(forecasts, axis)
    if forecasts.shape == observations.shape:
        forecasts = forecasts[..., np.newaxis]
    if observations.shape != forecasts.shape[:-1]:
        raise ValueError('observations and forecasts must have matching '
                         'shapes or matching shapes except along `axis=%s`'
                         % axis)
    scalar_threshold = threshold.ndim == 0
    if threshold.ndim > 1:
        raise ValueError('threshold must be scalar or 1-dimensional')
    if threshold.ndim == 1 and not (np.sort(threshold) == threshold).all():
        raise ValueError('1D thresholds must be sorted')
    threshold = threshold.reshape((1,) * observations.ndim + (-1,))
    if not issorted:
        forecasts = np.sort(forecasts, axis=-1)
    result = _threshold_brier_score_core(observations, forecasts, threshold)
    if scalar_threshold:
        result = result.squeeze(axis=-1)
    return result