0
# Statistical Downscaling and Bias Adjustment
1

2
Comprehensive statistical downscaling and bias correction methods for climate model post-processing. 
3

4
**Important Note**: The SDBA (Statistical Downscaling and Bias Adjustment) functionality has been moved to a separate `xsdba` package. While still accessible through `xclim.sdba` for backward compatibility, users are encouraged to use `xsdba` directly.
5

6
## Import Methods
7

8
**Legacy import (with deprecation warning):**
9
```python
10
from xclim import sdba
11
```
12

13
**Recommended import:**
14
```python
15
import xsdba as sdba
16
# or
17
from xsdba import EmpiricalQuantileMapping, DetrendedQuantileMapping
18
```
19

20
## Capabilities
21

22
### Bias Adjustment Methods
23

24
Core bias correction algorithms for correcting systematic errors in climate model simulations.
25

26
```python { .api }
27
class EmpiricalQuantileMapping:
28
    """
29
    Empirical Quantile Mapping bias adjustment method.
30
    
31
    Corrects the distribution of climate model data to match observed
32
    distribution using empirical quantile mapping.
33
    """
34
    
35
    def __init__(self, nquantiles=50, kind="additive", **kwargs):
36
        """
37
        Parameters:
38
        - nquantiles: int, number of quantiles for mapping (default 50)
39
        - kind: str, adjustment kind ("additive", "multiplicative")
40
        - **kwargs: additional method parameters
41
        """
42
    
43
    def train(self, ref, hist, **kwargs):
44
        """
45
        Train the bias adjustment model.
46
        
47
        Parameters:
48
        - ref: xr.DataArray, reference (observed) data for training
49
        - hist: xr.DataArray, historical model data for training
50
        - **kwargs: additional training parameters
51
        
52
        Returns:
53
        self: Trained bias adjustment object
54
        """
55
    
56
    def adjust(self, sim, **kwargs):
57
        """
58
        Apply bias adjustment to simulation data.
59
        
60
        Parameters:
61
        - sim: xr.DataArray, simulation data to adjust
62
        - **kwargs: additional adjustment parameters
63
        
64
        Returns:
65
        xr.DataArray: Bias-adjusted simulation data
66
        """
67

68
class DetrendedQuantileMapping:
69
    """
70
    Detrended Quantile Mapping for preserving climate change signals.
71
    
72
    Separates trend from variability, applies quantile mapping to 
73
    variability, then adds back the trend.
74
    """
75
    
76
    def __init__(self, nquantiles=50, kind="additive", **kwargs):
77
        """
78
        Parameters:
79
        - nquantiles: int, number of quantiles for mapping
80
        - kind: str, adjustment kind ("additive", "multiplicative")
81
        - **kwargs: additional method parameters
82
        """
83

84
class QuantileDeltaMapping:
85
    """
86
    Quantile Delta Mapping preserving relative changes.
87
    
88
    Adjusts model quantiles while preserving the relative change
89
    between historical and future periods.
90
    """
91
    
92
    def __init__(self, nquantiles=50, kind="multiplicative", **kwargs):
93
        """
94
        Parameters:
95
        - nquantiles: int, number of quantiles for mapping
96
        - kind: str, adjustment kind ("additive", "multiplicative") 
97
        - **kwargs: additional method parameters
98
        """
99

100
class PrincipalComponents:
101
    """
102
    Principal Component Analysis based bias adjustment.
103
    
104
    Uses PCA to adjust multivariate climate data while preserving
105
    spatial and temporal correlations.
106
    """
107
    
108
    def __init__(self, n_components=None, **kwargs):
109
        """
110
        Parameters:
111
        - n_components: int, number of principal components to retain
112
        - **kwargs: additional PCA parameters
113
        """
114

115
class LOCI:
116
    """
117
    Local Intensity Scaling bias adjustment method.
118
    
119
    Adjusts wet-day precipitation intensities based on local
120
    scaling factors.
121
    """
122
    
123
    def __init__(self, thresh="0.05 mm/day", **kwargs):
124
        """
125
        Parameters:
126
        - thresh: str, wet day threshold
127
        - **kwargs: additional LOCI parameters
128
        """
129

130
class Scaling:
131
    """
132
    Simple scaling bias adjustment method.
133
    
134
    Applies additive or multiplicative scaling factors
135
    based on mean bias between model and observations.
136
    """
137
    
138
    def __init__(self, kind="additive", **kwargs):
139
        """
140
        Parameters:
141
        - kind: str, scaling kind ("additive", "multiplicative")
142
        - **kwargs: additional scaling parameters
143
        """
144
```
145

146
### Performance Measures
147

148
Statistical measures for evaluating bias adjustment performance and model skill.
149

150
```python { .api }
151
def bias(obs, sim, **kwargs):
152
    """
153
    Calculate bias between observations and simulations.
154
    
155
    Parameters:
156
    - obs: xr.DataArray, observed data
157
    - sim: xr.DataArray, simulated data
158
    - **kwargs: additional calculation parameters
159
    
160
    Returns:
161
    xr.DataArray: Bias values (sim - obs)
162
    """
163

164
def relative_bias(obs, sim, **kwargs):
165
    """
166
    Calculate relative bias as percentage.
167
    
168
    Parameters:
169
    - obs: xr.DataArray, observed data
170
    - sim: xr.DataArray, simulated data
171
    - **kwargs: additional calculation parameters
172
    
173
    Returns:
174
    xr.DataArray: Relative bias values as percentage
175
    """
176

177
def correlation(obs, sim, **kwargs):
178
    """
179
    Calculate correlation coefficient between obs and sim.
180
    
181
    Parameters:
182
    - obs: xr.DataArray, observed data
183
    - sim: xr.DataArray, simulated data
184
    - **kwargs: additional calculation parameters
185
    
186
    Returns:
187
    xr.DataArray: Correlation coefficients
188
    """
189

190
def mae(obs, sim, **kwargs):
191
    """
192
    Calculate Mean Absolute Error.
193
    
194
    Parameters:
195
    - obs: xr.DataArray, observed data
196
    - sim: xr.DataArray, simulated data
197
    - **kwargs: additional calculation parameters
198
    
199
    Returns:
200
    xr.DataArray: Mean absolute error values
201
    """
202

203
def rmse(obs, sim, **kwargs):
204
    """
205
    Calculate Root Mean Square Error.
206
    
207
    Parameters:
208
    - obs: xr.DataArray, observed data
209
    - sim: xr.DataArray, simulated data
210
    - **kwargs: additional calculation parameters
211
    
212
    Returns:
213
    xr.DataArray: Root mean square error values
214
    """
215

216
def taylor_stats(obs, sim, **kwargs):
217
    """
218
    Calculate Taylor diagram statistics.
219
    
220
    Parameters:
221
    - obs: xr.DataArray, observed data
222
    - sim: xr.DataArray, simulated data
223
    - **kwargs: additional calculation parameters
224
    
225
    Returns:
226
    dict: Dictionary with correlation, standard deviation ratio, and RMSE
227
    """
228
```
229

230
### Detrending Methods
231

232
Trend removal and restoration methods for climate change signal preservation.
233

234
```python { .api }
235
class LoessDetrend:
236
    """
237
    LOESS (locally weighted regression) detrending method.
238
    
239
    Removes smooth trends using locally weighted polynomial regression.
240
    """
241
    
242
    def __init__(self, f=0.2, niter=1, **kwargs):
243
        """
244
        Parameters:
245
        - f: float, smoothing parameter (0-1)
246
        - niter: int, number of iterations
247
        - **kwargs: additional LOESS parameters
248
        """
249
    
250
    def fit(self, data, **kwargs):
251
        """
252
        Fit detrending model to data.
253
        
254
        Parameters:
255
        - data: xr.DataArray, input data with time dimension
256
        - **kwargs: additional fitting parameters
257
        
258
        Returns:
259
        self: Fitted detrending object
260
        """
261
    
262
    def detrend(self, data, **kwargs):
263
        """
264
        Remove trend from data.
265
        
266
        Parameters:
267
        - data: xr.DataArray, input data to detrend
268
        - **kwargs: additional detrending parameters
269
        
270
        Returns:
271
        xr.DataArray: Detrended data
272
        """
273
    
274
    def retrend(self, data, trend, **kwargs):
275
        """
276
        Add trend back to detrended data.
277
        
278
        Parameters:
279
        - data: xr.DataArray, detrended data
280
        - trend: xr.DataArray, trend to add back
281
        - **kwargs: additional parameters
282
        
283
        Returns:
284
        xr.DataArray: Retrended data
285
        """
286

287
class PolyDetrend:
288
    """
289
    Polynomial detrending method.
290
    
291
    Removes polynomial trends of specified degree.
292
    """
293
    
294
    def __init__(self, degree=1, **kwargs):
295
        """
296
        Parameters:
297
        - degree: int, polynomial degree for detrending
298
        - **kwargs: additional polynomial parameters
299
        """
300
```
301

302
### N-dimensional Utilities
303

304
Advanced bias adjustment utilities for multivariate and high-dimensional climate data.
305

306
```python { .api }
307
def stack_variables(datasets, rechunk=True):
308
    """
309
    Stack multiple climate variables for multivariate adjustment.
310
    
311
    Parameters:
312
    - datasets: dict, mapping variable names to xr.DataArray objects
313
    - rechunk: bool, whether to rechunk stacked data (default True)
314
    
315
    Returns:
316
    xr.DataArray: Stacked multivariate data array
317
    """
318

319
def unstack_variables(stacked_data, reference_datasets):
320
    """
321
    Unstack multivariate data back to individual variables.
322
    
323
    Parameters:
324
    - stacked_data: xr.DataArray, stacked multivariate data
325
    - reference_datasets: dict, reference datasets for unstacking template
326
    
327
    Returns:
328
    dict: Mapping variable names to unstacked xr.DataArray objects
329
    """
330

331
def apply_correction(x, correction_factors, **kwargs):
332
    """
333
    Apply correction factors to climate data.
334
    
335
    Parameters:
336
    - x: xr.DataArray, input data to correct
337
    - correction_factors: xr.DataArray, correction factors
338
    - **kwargs: additional application parameters
339
    
340
    Returns:
341
    xr.DataArray: Corrected data
342
    """
343

344
def interp_on_quantiles(newx, xq, yq, **kwargs):
345
    """
346
    Interpolate values based on quantile mapping.
347
    
348
    Parameters:
349
    - newx: array-like, new values to interpolate
350
    - xq: array-like, quantile values from source distribution
351
    - yq: array-like, quantile values from target distribution
352
    - **kwargs: interpolation parameters
353
    
354
    Returns:
355
    array-like: Interpolated values
356
    """
357
```
358

359
### Statistical Properties
360

361
Climate statistical property calculations for bias adjustment validation.
362

363
```python { .api }
364
def annual_cycle_amplitude(data, **kwargs):
365
    """
366
    Calculate amplitude of annual temperature cycle.
367
    
368
    Parameters:
369
    - data: xr.DataArray, daily or monthly climate data
370
    - **kwargs: additional calculation parameters
371
    
372
    Returns:
373
    xr.DataArray: Annual cycle amplitude
374
    """
375

376
def annual_cycle_phase(data, **kwargs):
377
    """
378
    Calculate phase of annual cycle (timing of maximum).
379
    
380
    Parameters:
381
    - data: xr.DataArray, daily or monthly climate data
382
    - **kwargs: additional calculation parameters
383
    
384
    Returns:
385
    xr.DataArray: Day of year of annual maximum
386
    """
387

388
def spell_length_distribution(data, threshold, op=">=", **kwargs):
389
    """
390
    Calculate distribution of spell lengths above/below threshold.
391
    
392
    Parameters:
393
    - data: xr.DataArray, daily climate data
394
    - threshold: float, threshold value for spell detection
395
    - op: str, comparison operator (">=", "<=", ">", "<")
396
    - **kwargs: additional parameters
397
    
398
    Returns:
399
    xr.DataArray: Spell length distribution statistics
400
    """
401

402
def transition_probabilities(data, threshold, **kwargs):
403
    """
404
    Calculate wet-to-dry and dry-to-wet transition probabilities.
405
    
406
    Parameters:
407
    - data: xr.DataArray, daily precipitation data
408
    - threshold: float, wet day threshold
409
    - **kwargs: additional calculation parameters
410
    
411
    Returns:
412
    dict: Dictionary with transition probability arrays
413
    """
414
```
415

416
## Usage Examples
417

418
### Basic Bias Adjustment
419

420
```python
421
import xarray as xr
422
from xclim import sdba
423

424
# Load climate data
425
ref = xr.open_dataset("observations.nc").pr  # Reference observations
426
hist = xr.open_dataset("model_historical.nc").pr  # Historical model
427
sim = xr.open_dataset("model_future.nc").pr  # Future projections
428

429
# Initialize bias adjustment method
430
QM = sdba.EmpiricalQuantileMapping(nquantiles=100, kind="multiplicative")
431

432
# Train the adjustment
433
QM.train(ref, hist)
434

435
# Apply adjustment to future projections
436
sim_adjusted = QM.adjust(sim)
437
```
438

439
### Advanced Detrended Quantile Mapping
440

441
```python
442
# Use detrended quantile mapping for temperature
443
tas_ref = xr.open_dataset("observations.nc").tas
444
tas_hist = xr.open_dataset("model_historical.nc").tas
445
tas_sim = xr.open_dataset("model_future.nc").tas
446

447
# Initialize detrended method
448
DQM = sdba.DetrendedQuantileMapping(
449
    nquantiles=50,
450
    kind="additive",
451
    trend_preservation="multiplicative"
452
)
453

454
# Train and apply
455
DQM.train(tas_ref, tas_hist)
456
tas_adjusted = DQM.adjust(tas_sim)
457
```
458

459
### Multivariate Bias Adjustment
460

461
```python
462
# Multivariate adjustment of temperature and precipitation
463
variables = {
464
    "tas": tas_hist,
465
    "pr": pr_hist
466
}
467
ref_vars = {
468
    "tas": tas_ref, 
469
    "pr": pr_ref
470
}
471

472
# Stack variables for multivariate processing
473
stacked_hist = sdba.stack_variables(variables)
474
stacked_ref = sdba.stack_variables(ref_vars)
475

476
# Apply PCA-based adjustment
477
PCA = sdba.PrincipalComponents(n_components=10)
478
PCA.train(stacked_ref, stacked_hist)
479

480
# Adjust future data
481
sim_vars = {"tas": tas_sim, "pr": pr_sim}
482
stacked_sim = sdba.stack_variables(sim_vars)
483
stacked_adjusted = PCA.adjust(stacked_sim)
484

485
# Unstack back to individual variables
486
adjusted_vars = sdba.unstack_variables(stacked_adjusted, sim_vars)
487
```
488

489
### Performance Evaluation
490

491
```python
492
from xclim.sdba import measures
493

494
# Calculate performance metrics
495
bias_val = measures.bias(tas_ref, tas_hist)
496
rmse_val = measures.rmse(tas_ref, tas_hist)  
497
corr_val = measures.correlation(tas_ref, tas_hist)
498

499
# Taylor diagram statistics
500
taylor_stats = measures.taylor_stats(tas_ref, tas_hist)
501
print(f"Correlation: {taylor_stats['correlation']}")
502
print(f"Std ratio: {taylor_stats['std_ratio']}")
503

504
# Before and after adjustment comparison
505
bias_before = measures.bias(tas_ref, tas_hist)
506
bias_after = measures.bias(tas_ref, tas_adjusted)
507
print(f"Bias reduction: {abs(bias_after) < abs(bias_before)}")
508
```
509

510
### Statistical Properties Analysis
511

512
```python
513
from xclim.sdba import properties
514

515
# Compare annual cycles
516
amp_obs = properties.annual_cycle_amplitude(tas_ref)
517
amp_model = properties.annual_cycle_amplitude(tas_adjusted)
518

519
# Check spell length distributions for precipitation
520
spell_dist_obs = properties.spell_length_distribution(
521
    pr_ref, threshold=1.0, op=">="
522
)
523
spell_dist_model = properties.spell_length_distribution(
524
    pr_adjusted, threshold=1.0, op=">="
525
)
526

527
# Transition probabilities for precipitation
528
trans_prob_obs = properties.transition_probabilities(pr_ref, threshold=1.0)
529
trans_prob_model = properties.transition_probabilities(pr_adjusted, threshold=1.0)
530
```