# Handling of missing data

deterministic probabilistic

Missing data must be flagged as “not a number” regardless of the utility used in the `evalhyd` stack. Missing data can be present both in the streamflow observations and/or the streamflow predictions. All time steps where either observations or predictions are missing are not taken into account in the computation of the metrics (a method sometimes referred to as pairwise deletion).

That is to say that these two calls to `evalhyd` yield the same result:

```>>> import evalhyd
>>> import numpy as np
>>> (
...     evalhyd.evald(np.array([7, 3, 3, np.nan, 5]), np.array([5, 4, 3, 5, np.nan]), ["KGE"])
...     == evalhyd.evald(np.array([7, 3, 3]), np.array([5, 4, 3]), ["KGE"])
... )
True
```
```> all(
+    evalhyd::evald(c(7, 3, 3, NA, 5), c(5, 4, 3, 5, NA), c("KGE"))[]
+    == evalhyd::evald(c(7, 3, 3), c(5, 4, 3), c("KGE"))[]
+ )
 TRUE
```

probabilistic-only

For a given site, `evalhyd` expects that the time indices in the streamflow predictions and in the streamflow observations correspond to the same datetime. This means that, when evaluating several lead times at once, the time series of streamflow observations will typically be longer than any of the time series of streamflow predictions. Indeed, missing predictions will exist for the longer lead times at the beginning of the time series, because the forecast datetime plus the lead time will not correspond to the first observed datetime (required for shorter lead times), but a later one. In turn, missing predictions will exist for the shorter lead times at the end of the time series, because the forecast datetime plus the lead time will not correspond to the last observed datetime (required for longer lead times), but an earlier one. For these datetimes, missing predictions must be filled with “not a number” value to make sure that the length of the observations and predictions are the same, so that a single time series of streamflow observations can be used across the prediction lead times evaluated.

To illustrate, let’s look at the simple example below where streamflow forecasts (only one ensemble member shown) are issued daily from the 1st until the 4th of January 2017 with a one-day, two-day, and three-day lead time.

```time index       0           1           2           3           4           5
datetime         2017-01-02  2017-01-03  2017-01-04  2017-01-05  2017-01-06  2017-01-07

observations     351         367         377         378         330         324

predictions
1-day lead time  312         335         358         342         NaN         NaN
2-day lead time  NaN         341         364         351         332         NaN
3-day lead time  NaN         NaN         361         358         327         327
```

For a complete forecast evaluation, the one-day lead time will require observations from the 2nd until the 5th of January, the two-day lead time from the 3rd until the 6th, and the three-day lead time from the 4th until the 7th. Since only one time series of streamflow observations is expected by `evalhyd`, it must start on the 2nd and end on the 7th of January, while the streamflow predictions must “flag” those time steps for which they are not making predictions as “not a number” (identified as `NaN` in the above example). Indeed, each forecast can only feature 4 predicted values if forecasts were issued only on 4 consecutive days.

This example translates into the following call to `evalhyd` (where the ensemble member shown above is duplicated four times for the only purpose of illustration):

```>>> import evalhyd
>>> import numpy as np
>>> res = evalhyd.evalp(
...     q_obs=np.array([[351, 367, 377, 378, 330, 324]]),
...     q_prd=np.array([[[[312, 335, 358, 342, np.nan, np.nan],
...                       [312, 335, 358, 342, np.nan, np.nan],
...                       [312, 335, 358, 342, np.nan, np.nan],
...                       [312, 335, 358, 342, np.nan, np.nan]],
...                      [[np.nan, 341, 364, 351, 332, np.nan],
...                       [np.nan, 341, 364, 351, 332, np.nan],
...                       [np.nan, 341, 364, 351, 332, np.nan],
...                       [np.nan, 341, 364, 351, 332, np.nan]],
...                      [[np.nan, np.nan, 361, 358, 327, 327],
...                       [np.nan, np.nan, 361, 358, 327, 327],
...                       [np.nan, np.nan, 361, 358, 327, 327],
...                       [np.nan, np.nan, 361, 358, 327, 327]]]]),
...     metrics=["CRPS_FROM_ECDF"]
... )
```
```> res <- evalhyd::evalp(
+     q_obs = rbind(c(351, 367, 377, 378, 330, 324)),
+     q_prd = aperm(
+         array(
+             data = c(
+                 rbind(c(312, 335, 358, 342, NA, NA),
+                       c(312, 335, 358, 342, NA, NA),
+                       c(312, 335, 358, 342, NA, NA),
+                       c(312, 335, 358, 342, NA, NA)),
+                 rbind(c(NA, 341, 364, 351, 332, NA),
+                       c(NA, 341, 364, 351, 332, NA),
+                       c(NA, 341, 364, 351, 332, NA),
+                       c(NA, 341, 364, 351, 332, NA)),
+                 rbind(c(NA, NA, 361, 358, 327, 327),
+                       c(NA, NA, 361, 358, 327, 327),
+                       c(NA, NA, 361, 358, 327, 327),
+                       c(NA, NA, 361, 358, 327, 327))
+             ),
+             dim = c(4, 6, 3, 1)
+         ),
+         perm = c(4, 3, 1, 2)
+     ),
+     metrics = c("CRPS_FROM_ECDF")
+ )
```