"""
Binary Forecast Verification
============================

This example demonstrates how to perform binary forecast verification.

If you are not familiar with the terms, you may wish to refer to the following
resources:
    
* `Verification of categorical predictands by Anna Ghelli of ECMWF <https://www.swpc.noaa.gov/sites/default/files/images/u30/Verification%20of%20Categorical%20Predictands.pdf>`_
* `Verification Measures by Weather Forecasting ... On-Line <https://www.wxonline.info/topics/verif2.html>`_
"""

###########################################################
# Definitions
# --------------------------------------------------------
#

########################################################################
# Import all required modules and methods:

# Python package to manage warning message
import warnings
# Python package for time calculations
import pandas as pd
# Python package for numerical calculations
import numpy as np
# Python package for xarrays to read and handle netcdf data
import xarray as xr
# swirlspy traditional binary verification package
from swirlspy.ver.crosstab import contingency
# swirlspy verification metric package
import swirlspy.ver.metric as mt

warnings.filterwarnings("ignore")

# Generating sample forecast and observed data
# Dimensions: x, y, time (3 x 3 x 4)
forecast_data = [
    [[7, 4, 6, 5], [4, 7, 9, 2], [3, 1, 6, 5]],
    [[4, 6, 8, 1], [11, 0, 3, 6], [0, 3, 6, 8]],
    [[5, 3, 7, 5], [7, 7, 6, 2], [3, 8, 9, 2]]
]

observed_data = [
    [[8, 7, 11, 4], [6, 7, 6, 0], [3, 3, 7, 3]],
    [[8, 8, 9, 5], [12, 1, 2, 8], [1, 4, 3, 12]],
    [[5, 6, 6, 2], [3, 6, 4, 7], [5, 7, 8, 2]]
]

# Creating xarrays
x = np.arange(3)
y = np.arange(3)
time = pd.date_range('1/1/2011', periods=4, freq='D')
forecast = xr.DataArray(
    forecast_data,
    coords=[('x', x), ('y', y), ('time', time)])
observed = xr.DataArray(
    observed_data,
    coords=[('x', x), ('y', y), ('time', time)])

# Applying the contingency function
cont = contingency(5, forecast, observed)

# Generating skill metrics
pod = mt.pod(cont)
far = mt.far(cont)
csi = mt.csi(cont)
freq_bias = mt.freq_bias(cont)
accuracy = mt.accuracy(cont)
pofd = mt.pofd(cont)
hss = mt.hss(cont)
ets = mt.ets(cont)
f1 = mt.f1Score(forecast, observed, average='weighted')

# Displaying results as pandas series
d = {
    'Hits': cont[0],
    'Misses': cont[1],
    'False Alarm': cont[2],
    'Correct Negative': cont[3],
    'Probability of Detection': pod,
    'False Alarm Ratio': far,
    'Critical Success Index': csi,
    'Frequency Bias': freq_bias,
    'Accuracy': accuracy,
    'Probability of False Detection': pofd,
    'Heidke Skill Score': hss,
    'Equitable Threat Score': ets,
    'f1 Score': f1
}

pd_series_d = pd.Series(d)
print(pd_series_d)