Skip to content

Surfaces

Surface grid computation and interpolation: builds the 3-D concentration surface array and provides vectorised bilinear interpolation for daily estimation.

wrtds.surfaces

Surface grid estimation and bilinear interpolation for WRTDS.

compute_surface_index(sample)

Compute the surface grid parameters from sample data.

The grid spans the range of observed log-discharge and decimal year, with fixed resolution: 14 log-Q levels and time steps of 1/16 year (~23 days).

Parameters:

Name Type Description Default
sample

Populated sample DataFrame (must have LogQ and DecYear).

 required

Returns:

Type Description

Dict with keys: bottom_logq, top_logq, step_logq,

n_logq, bottom_year, top_year, step_year,

n_year, logq_grid, year_grid.
Source code in wrtds/surfaces.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def compute_surface_index(sample):
    """Compute the surface grid parameters from sample data.

    The grid spans the range of observed log-discharge and decimal year,
    with fixed resolution: 14 log-Q levels and time steps of 1/16 year
    (~23 days).

    Args:
        sample: Populated sample DataFrame (must have ``LogQ`` and ``DecYear``).

    Returns:
        Dict with keys: ``bottom_logq``, ``top_logq``, ``step_logq``,
        ``n_logq``, ``bottom_year``, ``top_year``, ``step_year``,
        ``n_year``, ``logq_grid``, ``year_grid``.
    """
    bottom_logq = sample['LogQ'].min() - 0.05
    top_logq = sample['LogQ'].max() + 0.05
    n_logq = 14
    step_logq = (top_logq - bottom_logq) / (n_logq - 1)

    bottom_year = np.floor(sample['DecYear'].min())
    top_year = np.ceil(sample['DecYear'].max())
    step_year = 1.0 / 16.0
    n_year = round((top_year - bottom_year) / step_year) + 1

    return {
        'bottom_logq': bottom_logq,
        'top_logq': top_logq,
        'step_logq': step_logq,
        'n_logq': n_logq,
        'bottom_year': bottom_year,
        'top_year': top_year,
        'step_year': step_year,
        'n_year': n_year,
        'logq_grid': np.linspace(bottom_logq, top_logq, n_logq),
        'year_grid': np.linspace(bottom_year, top_year, n_year),
    }

estimate_surfaces(sample, surface_index, window_y=7.0, window_q=2.0, window_s=0.5, min_num_obs=100, min_num_uncen=50, edge_adjust=True)

Fit the censored regression at every grid point to build the surfaces array.

This is the main computational bottleneck of WRTDS: one MLE solve per grid point, typically ~7 000 total (14 log-Q levels x ~500 time steps).

Parameters:

Name Type Description Default
sample

Populated sample DataFrame with DecYear, LogQ, SinDY, CosDY, ConcLow, ConcHigh, Uncen.

 required
surface_index

Grid parameters from :func:compute_surface_index.

 required
window_y

Time half-window in years.

 7.0
window_q

Discharge half-window in log units.

 2.0
window_s

Season half-window in fraction of year.

 0.5
min_num_obs

Minimum observations with nonzero weight.

 100
min_num_uncen

Minimum uncensored observations with nonzero weight.

 50
edge_adjust

Expand time window near record edges.

 True

Returns:

Type Description

3-D numpy array of shape (n_logq, n_year, 3) where layer 0 is

yHat, layer 1 is SE, and layer 2 is ConcHat.
Source code in wrtds/surfaces.py 
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def estimate_surfaces(
    sample,
    surface_index,
    window_y=7.0,
    window_q=2.0,
    window_s=0.5,
    min_num_obs=100,
    min_num_uncen=50,
    edge_adjust=True,
):
    """Fit the censored regression at every grid point to build the surfaces array.

    This is the main computational bottleneck of WRTDS: one MLE solve per
    grid point, typically ~7 000 total (14 log-Q levels x ~500 time steps).

    Args:
        sample: Populated sample DataFrame with ``DecYear``, ``LogQ``,
            ``SinDY``, ``CosDY``, ``ConcLow``, ``ConcHigh``, ``Uncen``.
        surface_index: Grid parameters from :func:`compute_surface_index`.
        window_y: Time half-window in years.
        window_q: Discharge half-window in log units.
        window_s: Season half-window in fraction of year.
        min_num_obs: Minimum observations with nonzero weight.
        min_num_uncen: Minimum uncensored observations with nonzero weight.
        edge_adjust: Expand time window near record edges.

    Returns:
        3-D numpy array of shape ``(n_logq, n_year, 3)`` where layer 0 is
        ``yHat``, layer 1 is ``SE``, and layer 2 is ``ConcHat``.
    """
    logq_grid = surface_index['logq_grid']
    year_grid = surface_index['year_grid']
    n_logq = len(logq_grid)
    n_year = len(year_grid)

    # Pre-extract sample arrays for efficiency
    sample_data = {
        'DecYear': sample['DecYear'].values,
        'LogQ': sample['LogQ'].values,
        'SinDY': sample['SinDY'].values,
        'CosDY': sample['CosDY'].values,
        'ConcLow': sample['ConcLow'].values,
        'ConcHigh': sample['ConcHigh'].values,
        'Uncen': sample['Uncen'].values.astype(bool),
    }
    dec_year_arr = sample_data['DecYear']
    logq_arr = sample_data['LogQ']
    uncen_arr = sample_data['Uncen']
    record_start = float(dec_year_arr.min())
    record_end = float(dec_year_arr.max())

    surfaces = np.empty((n_logq, n_year, 3))

    for j in range(n_logq):
        for k in range(n_year):
            target_logq = logq_grid[j]
            target_year = year_grid[k]

            w = compute_weights(
                dec_year_arr, logq_arr, uncen_arr,
                target_year, target_logq,
                window_y, window_q, window_s,
                min_num_obs=min_num_obs,
                min_num_uncen=min_num_uncen,
                edge_adjust=edge_adjust,
                record_start=record_start,
                record_end=record_end,
            )

            beta, sigma = run_surv_reg(sample_data, w)
            yHat, SE, ConcHat = predict(beta, sigma, target_year, target_logq)

            surfaces[j, k, 0] = float(yHat)
            surfaces[j, k, 1] = float(SE)
            surfaces[j, k, 2] = float(ConcHat)

    return surfaces

make_interpolator(surfaces, surface_index, layer=2)

Create a RegularGridInterpolator for one surface layer.

Parameters:

Name Type Description Default
surfaces

3-D array (n_logq, n_year, 3).

 required
surface_index

Grid parameters dict.

 required
layer

Which layer to interpolate (0=yHat, 1=SE, 2=ConcHat).

 2

Returns:

Name Type Description
A

class:scipy.interpolate.RegularGridInterpolator instance.

Out-of-range queries are clamped to the nearest grid boundary.
Source code in wrtds/surfaces.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def make_interpolator(surfaces, surface_index, layer=2):
    """Create a ``RegularGridInterpolator`` for one surface layer.

    Args:
        surfaces: 3-D array ``(n_logq, n_year, 3)``.
        surface_index: Grid parameters dict.
        layer: Which layer to interpolate (0=yHat, 1=SE, 2=ConcHat).

    Returns:
        A :class:`scipy.interpolate.RegularGridInterpolator` instance.
        Out-of-range queries are clamped to the nearest grid boundary.
    """
    logq_grid = surface_index['logq_grid']
    year_grid = surface_index['year_grid']

    return RegularGridInterpolator(
        (logq_grid, year_grid),
        surfaces[:, :, layer],
        method='linear',
        bounds_error=False,
        fill_value=None,  # nearest-neighbour extrapolation
    )

interpolate_surface(surfaces, surface_index, logq, dec_year, layer=2)

Vectorised bilinear interpolation on a surface layer.

Queries are clamped to the grid boundaries before interpolation.

Parameters:

Name Type Description Default
surfaces

3-D array (n_logq, n_year, 3).

 required
surface_index

Grid parameters dict.

 required
logq

Log-discharge value(s).

 required
dec_year

Decimal year value(s).

 required
layer

Which layer (0=yHat, 1=SE, 2=ConcHat).

 2

Returns:

Type Description

Array of interpolated values (same shape as inputs).
Source code in wrtds/surfaces.py 
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def interpolate_surface(surfaces, surface_index, logq, dec_year, layer=2):
    """Vectorised bilinear interpolation on a surface layer.

    Queries are clamped to the grid boundaries before interpolation.

    Args:
        surfaces: 3-D array ``(n_logq, n_year, 3)``.
        surface_index: Grid parameters dict.
        logq: Log-discharge value(s).
        dec_year: Decimal year value(s).
        layer: Which layer (0=yHat, 1=SE, 2=ConcHat).

    Returns:
        Array of interpolated values (same shape as inputs).
    """
    logq = np.asarray(logq, dtype=float)
    dec_year = np.asarray(dec_year, dtype=float)

    # Clamp to grid boundaries
    logq_c = np.clip(logq, surface_index['logq_grid'][0], surface_index['logq_grid'][-1])
    dec_year_c = np.clip(dec_year, surface_index['year_grid'][0], surface_index['year_grid'][-1])

    interp = make_interpolator(surfaces, surface_index, layer)
    pts = np.column_stack([logq_c.ravel(), dec_year_c.ravel()])
    result = interp(pts)

    return result.reshape(logq.shape) if logq.ndim > 0 else result.item()