0
# Spatial Scores
1

2
Spatial verification methods that account for spatial structure and displacement errors in gridded forecasts. These metrics evaluate forecast performance while considering the spatial context and neighborhood relationships, addressing the "double penalty" problem in traditional point-wise verification.
3

4
## Capabilities
5

6
### Fractions Skill Score (FSS)
7

8
The primary spatial verification method using a neighborhood approach to evaluate forecast accuracy over different spatial scales.
9

10
#### 2D Fractions Skill Score
11

12
Standard FSS implementation for 2D gridded data with configurable spatial windows and event thresholds.
13

14
```python { .api }
15
def fss_2d(
16
    fcst: xr.DataArray,
17
    obs: xr.DataArray,
18
    *,
19
    event_threshold: float,
20
    window_size: Tuple[int, int],
21
    spatial_dims: Tuple[str, str],
22
    zero_padding: bool = False,
23
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
24
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
25
    threshold_operator: Optional[Callable] = None,
26
    compute_method: FssComputeMethod = FssComputeMethod.NUMPY,
27
    dask: str = "forbidden",
28
) -> xr.DataArray:
29
    """
30
    Calculate 2D Fractions Skill Score for spatial verification.
31

32
    Args:
33
        fcst: Forecast data array (2D spatial + other dims)
34
        obs: Observation data array with matching spatial dimensions
35
        event_threshold: Threshold for binary event definition
36
        window_size: (height, width) of spatial verification window
37
        spatial_dims: Names of spatial dimensions (e.g., ("lat", "lon"))
38
        zero_padding: Add zeros around domain edges for window calculations
39
        reduce_dims: Dimensions to reduce
40
        preserve_dims: Dimensions to preserve
41
        threshold_operator: Function for threshold comparison (default: >=)
42
        compute_method: Computational approach (NUMPY, SCIPY, etc.)
43
        dask: Dask computation control ("forbidden", "allowed")
44

45
    Returns:
46
        FSS values [0, 1] where 1 = perfect skill, 0 = no skill
47

48
    Formula:
49
        FSS = 1 - MSE_n / MSE_n_ref
50

51
    Where:
52
        - MSE_n: Mean squared error of neighborhood fractions
53
        - MSE_n_ref: Reference MSE (worst possible forecast)
54

55
    Notes:
56
        - Addresses the "double penalty" problem
57
        - Window size controls spatial scale of verification
58
        - Larger windows generally yield higher FSS values
59
        - Event threshold defines binary events from continuous data
60
    """
61
```
62

63
#### Binary FSS
64

65
FSS calculation for pre-discretized binary fields, skipping the threshold conversion step.
66

67
```python { .api }
68
def fss_2d_binary(
69
    fcst: xr.DataArray,
70
    obs: xr.DataArray,
71
    *,
72
    window_size: Tuple[int, int],
73
    spatial_dims: Tuple[str, str],
74
    zero_padding: bool = False,
75
    reduce_dims: Optional[FlexibleDimensionTypes] = None,
76
    preserve_dims: Optional[FlexibleDimensionTypes] = None,
77
    compute_method: FssComputeMethod = FssComputeMethod.NUMPY,
78
    dask: str = "forbidden",
79
) -> xr.DataArray:
80
    """
81
    Calculate FSS for binary data fields.
82

83
    Args:
84
        fcst: Binary forecast data {0, 1}
85
        obs: Binary observation data {0, 1}
86
        window_size: (height, width) of spatial verification window
87
        spatial_dims: Names of spatial dimensions
88
        zero_padding: Add zeros around domain edges
89
        reduce_dims: Dimensions to reduce
90
        preserve_dims: Dimensions to preserve
91
        compute_method: Computational approach
92
        dask: Dask computation control
93

94
    Returns:
95
        FSS values [0, 1]
96

97
    Notes:
98
        - Input data must already be binary {0, 1}
99
        - More efficient when events are pre-defined
100
        - Same FSS formula as fss_2d but skips thresholding
101
    """
102
```
103

104
#### Single Field FSS
105

106
FSS implementation optimized for single 2D spatial fields without additional dimensions.
107

108
```python { .api }
109
def fss_2d_single_field(
110
    fcst: xr.DataArray,
111
    obs: xr.DataArray,
112
    *,
113
    event_threshold: float,
114
    window_size: Tuple[int, int],
115
    spatial_dims: Tuple[str, str],
116
    zero_padding: bool = False,
117
    threshold_operator: Optional[Callable] = None,
118
    compute_method: FssComputeMethod = FssComputeMethod.NUMPY,
119
) -> xr.DataArray:
120
    """
121
    Calculate FSS for single 2D spatial field.
122

123
    Args:
124
        fcst: 2D forecast field
125
        obs: 2D observation field
126
        event_threshold: Threshold for event definition
127
        window_size: (height, width) of verification window
128
        spatial_dims: Names of spatial dimensions
129
        zero_padding: Add zeros around domain edges
130
        threshold_operator: Comparison function for thresholding
131
        compute_method: Computational approach
132

133
    Returns:
134
        Single FSS value
135

136
    Notes:
137
        - Optimized for single spatial fields
138
        - No additional dimensions (time, ensemble, etc.)
139
        - Returns scalar FSS value
140
    """
141
```
142

143
### FSS Computation Methods
144

145
Enumeration of available computational approaches for FSS calculation.
146

147
```python { .api }
148
from enum import Enum
149

150
class FssComputeMethod(Enum):
151
    """
152
    Available computation methods for FSS calculation.
153

154
    Attributes:
155
        NUMPY: Pure NumPy implementation (default)
156
        SCIPY: SciPy-based convolution approach
157
        DASK: Dask-distributed computation (if available)
158
    """
159
    NUMPY = "numpy"
160
    SCIPY = "scipy"
161
    DASK = "dask"
162
```
163

164
## Usage Patterns
165

166
### Basic Spatial Verification
167

168
```python
169
from scores.spatial import fss_2d
170
import xarray as xr
171
import numpy as np
172

173
# Create sample 2D gridded data
174
forecast = xr.DataArray(
175
    np.random.exponential(2, (50, 40)),  # 50x40 grid
176
    dims=["y", "x"],
177
    coords={"y": range(50), "x": range(40)}
178
)
179
observation = xr.DataArray(
180
    np.random.exponential(2, (50, 40)),
181
    dims=["y", "x"],
182
    coords={"y": range(50), "x": range(40)}
183
)
184

185
# Calculate FSS for precipitation >= 5mm with 5x5 window
186
fss = fss_2d(
187
    forecast, observation,
188
    event_threshold=5.0,
189
    window_size=(5, 5),
190
    spatial_dims=("y", "x")
191
)
192

193
print(f"FSS (5x5 window, 5mm threshold): {fss.values:.3f}")
194
```
195

196
### Multi-scale Spatial Analysis
197

198
```python
199
# Evaluate FSS across multiple spatial scales
200
window_sizes = [(1, 1), (3, 3), (5, 5), (11, 11), (21, 21)]
201
fss_scales = {}
202

203
for window in window_sizes:
204
    fss_value = fss_2d(
205
        forecast, observation,
206
        event_threshold=5.0,
207
        window_size=window,
208
        spatial_dims=("y", "x")
209
    )
210
    scale_name = f"{window[0]}x{window[1]}"
211
    fss_scales[scale_name] = fss_value.values
212

213
print("FSS by spatial scale:")
214
for scale, fss_val in fss_scales.items():
215
    print(f"  {scale}: {fss_val:.3f}")
216
```
217

218
### Multi-threshold Analysis
219

220
```python
221
# Evaluate FSS for different event intensities
222
thresholds = [1.0, 2.5, 5.0, 10.0, 20.0]
223
window_size = (7, 7)
224

225
fss_thresholds = {}
226
for threshold in thresholds:
227
    fss_value = fss_2d(
228
        forecast, observation,
229
        event_threshold=threshold,
230
        window_size=window_size,
231
        spatial_dims=("y", "x")
232
    )
233
    fss_thresholds[threshold] = fss_value.values
234

235
print(f"FSS by threshold (window: {window_size[0]}x{window_size[1]}):")
236
for thresh, fss_val in fss_thresholds.items():
237
    print(f"  {thresh:4.1f}mm: {fss_val:.3f}")
238
```
239

240
### Time Series Spatial Verification
241

242
```python
243
# Multi-temporal spatial verification
244
forecast_time = xr.DataArray(
245
    np.random.exponential(2, (24, 50, 40)),  # 24 time steps
246
    dims=["time", "y", "x"],
247
    coords={
248
        "time": pd.date_range("2023-01-01", periods=24, freq="H"),
249
        "y": range(50),
250
        "x": range(40)
251
    }
252
)
253
observation_time = xr.DataArray(
254
    np.random.exponential(2, (24, 50, 40)),
255
    dims=["time", "y", "x"],
256
    coords=forecast_time.coords
257
)
258

259
# FSS for each time step
260
fss_hourly = fss_2d(
261
    forecast_time, observation_time,
262
    event_threshold=5.0,
263
    window_size=(7, 7),
264
    spatial_dims=("y", "x"),
265
    reduce_dims=["y", "x"]  # Reduce spatial dims, keep time
266
)
267

268
print(f"Hourly FSS shape: {fss_hourly.shape}")
269
print(f"Mean FSS over 24 hours: {fss_hourly.mean().values:.3f}")
270

271
# Overall FSS (all time steps combined)
272
fss_overall = fss_2d(
273
    forecast_time, observation_time,
274
    event_threshold=5.0,
275
    window_size=(7, 7),
276
    spatial_dims=("y", "x")
277
    # No reduce_dims specified = reduce all non-spatial dims
278
)
279

280
print(f"Overall FSS: {fss_overall.values:.3f}")
281
```
282

283
### Binary Data FSS
284

285
```python
286
from scores.processing import binary_discretise
287

288
# Pre-process continuous data to binary events
289
precip_threshold = 10.0
290
fcst_binary = binary_discretise(forecast_time, precip_threshold)
291
obs_binary = binary_discretise(observation_time, precip_threshold)
292

293
# Calculate FSS for binary data (more efficient)
294
fss_binary = fss_2d_binary(
295
    fcst_binary, obs_binary,
296
    window_size=(9, 9),
297
    spatial_dims=("y", "x")
298
)
299

300
print(f"FSS for binary data (10mm threshold): {fss_binary.values:.3f}")
301
```
302

303
### Advanced FSS Options
304

305
```python
306
import operator
307

308
# Custom threshold operator (less than or equal)
309
fss_le = fss_2d(
310
    forecast, observation,
311
    event_threshold=2.0,
312
    window_size=(5, 5),
313
    spatial_dims=("y", "x"),
314
    threshold_operator=operator.le  # Event = value <= 2.0
315
)
316

317
# With zero padding (extends domain with zeros)
318
fss_padded = fss_2d(
319
    forecast, observation,
320
    event_threshold=5.0,
321
    window_size=(11, 11),
322
    spatial_dims=("y", "x"),
323
    zero_padding=True  # Adds zeros around edges
324
)
325

326
# Different computation method
327
fss_scipy = fss_2d(
328
    forecast, observation,
329
    event_threshold=5.0,
330
    window_size=(7, 7),
331
    spatial_dims=("y", "x"),
332
    compute_method=FssComputeMethod.SCIPY
333
)
334

335
print(f"FSS (<=2.0mm events): {fss_le.values:.3f}")
336
print(f"FSS (with padding): {fss_padded.values:.3f}")
337
print(f"FSS (SciPy method): {fss_scipy.values:.3f}")
338
```
339

340
### Performance Considerations
341

342
```python
343
# For large datasets, control computation methods
344
large_forecast = xr.DataArray(
345
    np.random.exponential(2, (100, 200, 300)),  # Large 3D array
346
    dims=["time", "y", "x"]
347
)
348
large_observation = xr.DataArray(
349
    np.random.exponential(2, (100, 200, 300)),
350
    dims=["time", "y", "x"]
351
)
352

353
# Process in chunks for memory efficiency
354
fss_chunked = fss_2d(
355
    large_forecast.chunk({"time": 10}),  # Chunk along time
356
    large_observation.chunk({"time": 10}),
357
    event_threshold=5.0,
358
    window_size=(7, 7),
359
    spatial_dims=("y", "x"),
360
    reduce_dims=["y", "x"],  # Keep time dimension
361
    compute_method=FssComputeMethod.NUMPY,
362
    dask="allowed"  # Allow Dask computation
363
)
364

365
# Compute results
366
fss_result = fss_chunked.compute()  # Triggers Dask computation
367
print(f"FSS time series computed: shape {fss_result.shape}")
368
```
369

370
## Interpretation Guidelines
371

372
### FSS Values and Skill
373

374
- **FSS = 1.0**: Perfect spatial forecast
375
- **FSS = 0.5**: Useful skill threshold (commonly used benchmark)
376
- **FSS = 0.0**: No skill (random forecast)
377
- **FSS < 0.0**: Forecast worse than climatology (rare for reasonable forecasts)
378

379
### Window Size Effects
380

381
- **Small windows** (1x1 to 3x3): Strict point-wise verification
382
- **Medium windows** (5x5 to 11x11): Moderate spatial tolerance
383
- **Large windows** (21x21+): Very permissive spatial verification
384

385
### Threshold Sensitivity
386

387
- **Low thresholds**: Easy events, higher FSS values
388
- **High thresholds**: Rare events, lower FSS values, more challenging verification
389

390
### Scale-dependent Interpretation
391

392
FSS naturally increases with window size, so comparisons should use consistent spatial scales or reference the "useful scale" where FSS first exceeds 0.5.