Point-to-Grid Matching

Comparing a model against point observations at the model's own nearest cell punishes a forecast that is displaced just as hard as one that is wrong. Relaxing the comparison to "the best-matching cell within R" separates those two failures — the model got the amount right but put it in the wrong place, versus the model got the amount wrong.

A neighborhood search is a biased estimator, always

In a 3 km box on a 1 km grid there are roughly 50 candidate cells. A near-matching value exists almost regardless of whether the model has any skill at that location, so best-matching always improves the apparent fit, including for a model with no skill at all.

An improvement figure quoted on its own is therefore not a result. It has to be compared against null_improvement — the same search re-run with the whole point set displaced to deliberately wrong locations. If the real improvement sits inside that distribution, what you measured was the search's optimism, not the model.

Which function answers which question

function	question
neighborhood_match	Could the model have been right nearby? How intense does it get nearby at all?
best_match_locate	Where is the match, and when did it fall there?
grid_best_match	Values only, for a whole point set — the null's workhorse
logvar_improvement	How much scatter did the search remove?
null_improvement	How much of that was luck?

best_match_locate is the one that turns matching into a measurement: knowing a matching value sits 4 km away is weaker than knowing the match sits consistently northwest of every gauge. The first is a scatter; the second is a displacement vector.

Why the null uses one common offset

null_improvement displaces every point by the same random offset each trial, rather than jittering each point independently. That preserves the network's geometry and its relationship to the field's spatial structure, and destroys only the correspondence — which is what makes it a null for displacement specifically.

Do not restructure the null's sampling loop

The generator is consumed inside a rejection loop, so the number of draws depends on the data. Each attempt takes exactly a bearing then a distance, including attempts later discarded for landing off-footprint. Batching the draws with a size= argument, reordering them, or precomputing the offsets changes every downstream number while looking like a speed-up.

Example

import numpy as np
from modverif.match import grid_best_match, logvar_improvement, null_improvement

radius_m = 3000.0
point_vals, best_vals = grid_best_match(field, gx, gy, sx, sy, obs, radius_m)
real, n_used = logvar_improvement(point_vals, best_vals, obs)

null = null_improvement(field, gx, gy, sx, sy, obs, radius_m,
                        n_trials=999, rng=np.random.default_rng(0))

if len(null):                       # may be empty -- percentile of an empty array raises
    p = (np.sum(null >= real) + 1) / (len(null) + 1)
    print(f'improvement {real:.3f} vs null median {np.median(null):.3f}, p = {p:.3f}')

Coordinates must be projected

All coordinates and radii are in a projected, metric CRS — the boxes are built by searchsorted on each axis in the axis's own units. Passing degrees produces a box whose east–west extent varies with latitude.

Boxes are square (Chebyshev), not circular. That is deliberate: it makes each search a pair of cheap axis slices rather than a distance computation over every candidate cell.

API

Point-to-grid neighborhood matching, and the selection-bias null that keeps it honest.

When a model is compared against point observations at the model's own nearest cell, a spatially displaced but otherwise correct forecast scores badly. Relaxing the comparison to "the best-matching cell within R" separates an amplitude error from a placement error -- but only if you also price what that search buys by luck, which is what null_improvement is for.

NOTE: This is the point-based neighborhood family. The grid-based one lives in modverif.metrics -- compute_fss, compute_fraction_field, compute_fss_multi_scale. The fractions skill score compares two fields by the fraction of each neighborhood exceeding a threshold; the functions here match scattered points against a field by value. Related ideas, different inputs, different answers.

WARNING: Searching a neighborhood for the closest value is a biased estimator, always. In a 3 km box on a 1 km grid there are ~50 candidate cells, so a near-matching value exists almost regardless of whether the model has any real skill at that location. Any improvement from best-matching must be compared against null_improvement -- the same search re-run at deliberately wrong locations -- or it means nothing. Reporting the improvement alone overstates the model.

Coordinates are in a projected, metric CRS; radii are in the same units (metres). Boxes are square (Chebyshev), not circular, and that is deliberate: it makes the search a cheap searchsorted slice on each axis.

best_match_locate(field, start_hours, gx, gy, px, py, obs, radius_m)

Locate the best-matching cell for a single point: its value, where it is, and its timing.

The extra return values are the point. Knowing that a matching value exists 4 km away is a weaker statement than knowing the match sits consistently to the northwest of every gauge, which is what turns a scatter of best-matches into a displacement measurement.

Parameters:

Name	Type	Description	Default
field	ndarray	Gridded model values, shape (len(gy), len(gx)).	required
start_hours	ndarray	Per-cell accumulation-window start, same shape as field -- so the matched cell's timing comes for free rather than needing a second lookup.	required
gx	ndarray	Ascending grid axis coordinates, projected metric CRS.	required
gy	ndarray	Ascending grid axis coordinates, projected metric CRS.	required
px	float	Point coordinates in the same CRS.	required
py	float	Point coordinates in the same CRS.	required
obs	float	Observed value to match.	required
radius_m	float	Half-width of the square search box.	required

Returns:

Name	Type	Description
value	float	Best-matching cell value, NaN if the box held no finite cell.
	iy, ix : int	Indices of that cell in field, or -1 if none.
start	float	That cell's window start, NaN if none.

field_shift_objective(field, gx, gy, sx, sy, obs, max_shift_km, step_km)

One whole-field displacement vector, with the objective surface that produced it.

Where neighborhood_match lets every point find its own best cell -- and can therefore cherry-pick -- this asks a single question of the entire field: which rigid translation makes the model best match all the points at once? One number, no per-point freedom, and the surface around the optimum shows whether it is well determined or a broad flat basin.

The scored point set is fixed across every candidate shift: only points whose full ±max_shift box is finite and strictly positive take part. Otherwise shifts that happen to move points onto the footprint would score against a different sample from those that do not, and the comparison would be meaningless.

Parameters:

Name	Type	Description	Default
field	ndarray	Gridded model values, shape (len(gy), len(gx)).	required
gx	ndarray	Ascending, regularly spaced grid axes in a projected metric CRS.	required
gy	ndarray	Ascending, regularly spaced grid axes in a projected metric CRS.	required
sx	ndarray	Point coordinates in the same CRS.	required
sy	ndarray	Point coordinates in the same CRS.	required
obs	ndarray	Observed value per point; only strictly positive values participate (the objective is a log-ratio).	required
max_shift_km	float	Half-width of the shift search, km.	required
step_km	float	Shift increment, km; rounded to at least one grid cell.	required

Returns:

Name	Type	Description
offsets_km	ndarray or None	Shift values along each axis, km. None if fewer than MIN_GAUGES_FOR_IMPROVEMENT points survived the fixed-sample rule.
objective	ndarray or None	J[a, b] = mean squared log-ratio at shift (dx=offsets[b], dy=offsets[a]) -- rows are northing, columns are easting. None when offsets_km is None.
used	ndarray	Boolean mask of the points that were scored, returned even on the None path so a caller can report why it declined.

Notes

Sign convention. A positive optimum (u, v) means the model values that match the observations sit at point + (u, v) -- i.e. the field is displaced by +(u, v). Correcting it means shifting the field by -(u, v). Getting this backwards inverts every conclusion drawn from it, so it is stated here and in the returned array's axis order.

grid_best_match(field, gx, gy, sx, sy, obs, radius_m)

Point values and best-match values for a whole point set at one radius.

The values-only workhorse behind null_improvement, which calls it once per null trial.

Parameters:

Name	Type	Description	Default
field	ndarray	Gridded model values, shape (len(gy), len(gx)).	required
gx	ndarray	Ascending, regularly spaced grid axes (see nearest_indices).	required
gy	ndarray	Ascending, regularly spaced grid axes (see nearest_indices).	required
sx	ndarray	Point coordinates in the same CRS.	required
sy	ndarray	Point coordinates in the same CRS.	required
obs	ndarray	Observed value per point.	required
radius_m	float	Half-width of the square search box.	required

Returns:

Name	Type	Description
point_vals	ndarray	Value at each point's own nearest cell, NaN where that cell is not finite.
best_vals	ndarray	Closest value to obs within the box, NaN where the box held no finite cell.

logvar_improvement(point_vals, best_vals, obs)

Fractional reduction in log-ratio variance from allowing the neighborhood search.

1 - var(log(best / obs)) / var(log(point / obs)).

Fractional, not absolute, and that is what makes the null comparison fair. Random locations start from a much worse point fit, so an absolute difference would flatter them mechanically; a ratio asks the same question of the real and displaced gauge sets.

Parameters:

Name	Type	Description	Default
point_vals	ndarray	Nearest-cell and best-match values, as from grid_best_match.	required
best_vals	ndarray	Nearest-cell and best-match values, as from grid_best_match.	required
obs	ndarray	Observed value per point.	required

Returns:

Name	Type	Description
improvement	float	Fraction of log-variance removed, or NaN if fewer than MIN_GAUGES_FOR_IMPROVEMENT points are usable or the point fit had no variance.
n_used	int	How many points were finite and strictly positive in all three inputs -- the log needs positives, so zeros and gaps drop out.

nearest_indices(coords, points)

Nearest index per point on an ascending, regularly spaced axis.

Computed arithmetically from the first two coordinates rather than by search, so it is O(1) per point -- and therefore assumes even spacing. On an irregular axis use np.searchsorted.

Parameters:

Name	Type	Description	Default
coords	ndarray	1-D ascending, regularly spaced axis.	required
points	ndarray	Positions to locate, same units.	required

Returns:

Type	Description
ndarray	Integer indices, clipped to the axis bounds -- a point outside the axis snaps to the nearest end rather than raising.

neighborhood_match(field, gx, gy, sx, sy, obs, radii_m, within_frac=0.2)

Per point and per search radius: the best-matching cell value, the neighborhood maximum, and whether any cell falls within a tolerance of the observation.

The three answer different questions. bm asks "could the model have been right nearby?", nmx asks "how intense does the model get nearby at all?", and near reduces the first to a pass/fail that survives tabulation.

Parameters:

Name	Type	Description	Default
field	ndarray	Gridded model values, shape (len(gy), len(gx)). NaN marks off-footprint cells.	required
gx	ndarray	Ascending grid axis coordinates in a projected, metric CRS.	required
gy	ndarray	Ascending grid axis coordinates in a projected, metric CRS.	required
sx	ndarray	Point coordinates in the same CRS.	required
sy	ndarray	Point coordinates in the same CRS.	required
obs	ndarray	Observed value per point.	required
radii_m	iterable of float	Half-widths of the square search boxes, in the CRS's units.	required
within_frac	float	Relative tolerance for the pass/fail flag.	0.2

Returns:

Type	Description
dict	{radius: {'bm': ndarray, 'nmx': ndarray, 'near': ndarray}}, each array aligned to the input points. Entries are NaN (or False) where the box contained no finite cell.

null_improvement(field, gx, gy, sx, sy, obs, radius_m, n_trials, rng, min_frac=0.8, null_dist_km=NULL_DIST_KM)

Price what the neighborhood search buys by luck alone.

Per trial, displace every point by one common random offset and recompute the point-to-best-match improvement there. A common offset rather than independent per-point jitter is what makes this a null for displacement: it preserves the gauge network's geometry and its relationship to the field's spatial structure, and destroys only the correspondence.

Compare the real improvement against this distribution. If it sits inside, the apparent skill of best-matching is search luck.

WARNING: The generator is consumed inside a rejection loop, so the number of draws depends on the data. Each attempt takes exactly two values (bearing, then distance) including attempts later discarded for landing off-footprint. Batching the draws, reordering them, or precomputing offsets changes every downstream number while looking like a speed-up. Pass a generator in; do not restructure this loop.

Parameters:

Name	Type	Description	Default
field	ndarray	Gridded model values.	required
gx	ndarray	Ascending, regularly spaced grid axes.	required
gy	ndarray	Ascending, regularly spaced grid axes.	required
sx	ndarray	True point coordinates -- displaced internally, not modified.	required
sy	ndarray	True point coordinates -- displaced internally, not modified.	required
obs	ndarray	Observed value per point.	required
radius_m	float	Half-width of the square search box; must match the radius being tested.	required
n_trials	int	Target number of usable trials. The loop gives up after 10 * n_trials attempts, so a short return means the footprint was hard to land on.	required
rng	Generator	Caller's generator. See the warning above.	required
min_frac	float	Minimum fraction of displaced points that must land on finite cells for a trial to count.	0.8
null_dist_km	tuple of float	(min, max) offset magnitude, km, sampled uniformly.	NULL_DIST_KM

Returns:

Type	Description
ndarray	Improvement per usable trial. May be shorter than n_trials, and may be empty -- take percentiles only after checking, since np.percentile of an empty array raises.

vector_coherence(dx, dy)

Do the offset vectors point the same way? The second selection-bias guard.

A real displacement makes neighbouring points' offsets align; search luck scatters them uniformly. This reduces that to one number -- the mean resultant length of the bearings, which is 1 for perfect alignment and near 0 for random directions -- plus a Rayleigh test of uniformity.

Use it alongside null_improvement. That one prices the magnitude of the apparent improvement; this one asks whether the offsets have a coherent direction, which a value-search on an unrelated field has no reason to produce.

Parameters:

Name	Type	Description	Default
dx	ndarray	Offset components per point, in a projected metric CRS. Pass only points that actually moved -- zero-length vectors have no bearing, and including them deflates the resultant.	required
dy	ndarray	Offset components per point, in a projected metric CRS. Pass only points that actually moved -- zero-length vectors have no bearing, and including them deflates the resultant.	required

Returns:

Name	Type	Description
rbar	float	Mean resultant length in [0, 1].
p	float	Rayleigh p-value against the uniform null, exp(-n * rbar**2).
	mean_dx, mean_dy : float	Mean offset vector, same units as the inputs.

Notes

Bearings are computed in compass convention -- arctan2(dx, dy), so 0 is north and angles increase clockwise. Note that rbar and p are invariant to that choice: rotating every bearing by a constant leaves the resultant length unchanged, so no convention is recoverable from them. The convention is documented because it governs how to interpret the returned mean offset vector, and because this is deliberately not the meteorological wind convention used by modverif.metrics.compute_wind_direction_bias -- the two answer different questions.