src.tests.raster_grid_postproc_ref

author:: Alberto M. Esmoris Pena

TEST-ONLY reference implementations of the five new raster post- processors introduced by the RasterGridEvaluatorPP task.

These pure-Python references act as cross-validation ground truth for the unit tests in src.tests.raster_grid_evaluatorpp_test.

IMPORTANT: production code must NEVER import from this module. The production C++/Python orchestrator lives in src/eval/raster/raster_postproc.py.

The implementations follow:

hillshade_horn_ref: Horn (1981) 3x3 derivative formulas in the pure GDAL math-CCW convention (plan section 4.2.1).
lrm_ref: Hesse (2010) local relief model.
resampling_filter_ref: SAGA/RVT downsample + bicubic spline upsample.
vat_ref: Kokalj and Somrak (2019) visualization for archaeological topography, RVT-style overlay + opacity blend.
msrm_ref: Orengo and Petrie (2018) multi-scale relief model, verbatim port of MSRM.ipynb’s compute_i_n + rm_task.

Functions

`hillshade_horn_ref`(z2d, xres, yres[, ...])	Compute the GDAL math-CCW Horn (1981) hillshade.
`lrm_ref`(z2d, radius_cells[, filter_kind, ...])	Compute the local relief model (Hesse, 2010).
`msrm_index_formula_ref`(xres, yres, fmin, fmax)	Compute the MSRM scale indices `(i, n)` and the radii list.
`msrm_ref`(z2d, xres, yres, fmin, fmax[, x, ...])	Multi-Scale Relief Model (Orengo & Petrie, 2018) verbatim port of `MSRM.ipynb`.
`resampling_filter_ref`(z2d, scale_factor)	Apply the SAGA/RVT resampling filter (downsample then bicubic-spline upsample).
`svf_po_cpp_match_ref`(z2d, xres, yres, ...)	Compute SVF and PO with the EXACT semantics of the C++ VAT kernel (`cpp/src/raster/RasterPostProcessors.tpp::vat`).
`svf_po_ref`(z2d, xres, yres, ...)	Compute sky-view factor and positive openness on a 2D DEM.
`uniform_filter_nan`(arr, size[, mode])	Compute a NaN-aware uniform (box-mean) filter over a 2D array.
`vat_ref`(z2d, xres, yres[, hillshade_az, ...])	Compute the RVT-style VAT composite (Kokalj & Somrak, 2019).

src.tests.raster_grid_postproc_ref.uniform_filter_nan(arr, size, mode='reflect')

Compute a NaN-aware uniform (box-mean) filter over a 2D array.

The numerator is the uniform filter of the NaN-zeroed input; the denominator is the uniform filter of the finite mask. The cell-wise mean is numerator / denominator (NaN where the denominator is zero).

This helper is part of the PUBLIC surface of this test-only reference module (no leading underscore) so it can be imported directly by the test suite.

Parameters:

arr (np.ndarray) – Input 2D array.
size (int) – Filter window full side length (must be odd).
mode (str) – Edge handling mode. Accepts the public "reflect" / "nearest" vocabulary (mapped to scipy’s flavour via _cpp_edge_to_scipy_mode()) so the reference bit-matches the C++ EdgeMode enum.

Returns:

The NaN-aware uniform-filtered array.

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.hillshade_horn_ref(z2d, xres, yres, azimuth_deg=315.0, altitude_deg=45.0, z_factor=1.0)

Compute the GDAL math-CCW Horn (1981) hillshade.

For each 3x3 patch [a b c; d e f; g h i] the Horn derivatives are:

\[ \begin{align}\begin{aligned}\frac{dz}{dx} = \frac{(c + 2f + i) - (a + 2d + g)}{8\,xres}\\\frac{dz}{dy} = \frac{(g + 2h + i) - (a + 2b + c)}{8\,yres}\end{aligned}\end{align} \]

The hillshade is then:

\[hs = \cos(\theta_z) \cos(\theta_s) + \sin(\theta_z) \sin(\theta_s) \cos(\theta_{az} - \theta_{asp})\]

clipped to \([0, 1]\). theta_{az} and theta_{asp} are both expressed in the same math-CCW convention so the formula matches gdaldem hillshade -alg Horn to fp32 noise.

Parameters:

z2d (np.ndarray) – 2D float DEM matrix.
xres (float) – Cell size along the x-axis.
yres (float) – Cell size along the y-axis.
azimuth_deg (float) – Solar azimuth in degrees clockwise from north (GDAL convention).
altitude_deg (float) – Solar altitude in degrees above the horizon.
z_factor (float) – Vertical exaggeration factor.

Returns:

Hillshade in [0, 1] with the same shape as z2d.

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.lrm_ref(z2d, radius_cells, filter_kind='mean', edge_mode='reflect')

Compute the local relief model (Hesse, 2010).

The LRM is z - lowpass(z, radius_cells) where the lowpass is either a box-mean (filter_kind='mean') or a Gaussian (filter_kind='gaussian', with sigma = radius / 3).

Parameters:

z2d (np.ndarray) – 2D DEM array.
radius_cells (int) – Lowpass radius in cells (window side is 2 * radius_cells + 1 for the mean filter).
filter_kind (str) – Either 'mean' or 'gaussian'.
edge_mode (str) – Boundary mode for the underlying scipy filter.

Returns:

The LRM map.

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.resampling_filter_ref(z2d, scale_factor)

Apply the SAGA/RVT resampling filter (downsample then bicubic-spline upsample).

Parameters:

z2d (np.ndarray) – 2D DEM array.
scale_factor (int) – Integer downsampling factor (must divide both dimensions).

Returns:

The resampled grid with the same shape as z2d.

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.svf_po_ref(z2d, xres, yres, search_radius_cells, n_directions)

Compute sky-view factor and positive openness on a 2D DEM.

Exposed on the public surface of the test-only reference module so subtests N07/N08 can call it directly.

Parameters:

z2d (np.ndarray) – 2D DEM array.
xres (float) – Cell size along x.
yres (float) – Cell size along y.
search_radius_cells (int) – Max ray walk length, in cells.
n_directions (int) – Number of azimuth rays per cell.

Returns:

Tuple (svf, po) of 2D arrays.

Return type:

tuple

src.tests.raster_grid_postproc_ref.svf_po_cpp_match_ref(z2d, xres, yres, search_radius_cells, n_directions)

Compute SVF and PO with the EXACT semantics of the C++ VAT kernel (cpp/src/raster/RasterPostProcessors.tpp::vat).

The pure-numerical companion svf_po_ref() implements the canonical RVT-style algorithm (max-elevation-angle horizon for SVF; in-bounds-skip boundary). The C++ kernel differs in two ways that matter at the boundary and on rough terrain:

Boundary: rays that step off the grid are reflect-clamped via the kernel’s EdgeMode::REFLECT (numpy np.pad(..., 'reflect') semantics, equivalent to scipy mode='mirror'), not skipped.
SVF math: the C++ stores alphaMax only on cells where v > zMax (i.e., the elevation angle taken at the highest-z cell along the ray, not the true max-elevation-angle cell). The two can differ; e.g., a far cell with a slightly higher elevation can win the highest-z tiebreak even though a closer cell carries a larger elevation angle.

Both choices are encoded as-is here so svf_po_ref() stays available for documentation / cross-reference uses (notebooks, RVT parity checks) while the test suite asserts against the C++-faithful version below.

Parameters:

z2d (np.ndarray) – 2D DEM array.
xres (float) – Cell size along x.
yres (float) – Cell size along y.
search_radius_cells (int) – Max ray walk length, in cells.
n_directions (int) – Number of azimuth rays per cell.

Returns:

Tuple (svf, po) of 2D arrays. The output policy matches the C++ kernel: PO is the un-normalised mean angle in radians, SVF is 1 - mean(sin(alphaMax)) clamped to [0, 1] by the kernel’s final blend.

Return type:

tuple

src.tests.raster_grid_postproc_ref.vat_ref(z2d, xres, yres, hillshade_az=315.0, hillshade_alt=45.0, hillshade_z=1.0, search_radius_cells=20, n_directions=16, svf_opacity=0.5, po_opacity=0.5)

Compute the RVT-style VAT composite (Kokalj & Somrak, 2019).

The composite blends a hillshade, a (1 - slope) layer, a sky-view factor layer and a positive-openness layer using multiplicative overlay + opacity blends.

Slope-unit convention. The slope layer is computed end-to-end in DEGREES on the Python side. The JSON schema and Python orchestrator carry slope thresholds in degrees; the C++ kernel performs its own deg -> rad conversion via the standard *pi/180 inside its internal kernel. Tests passing non-default values will get a one-time warning behaviour folded in by the production wrapper, not this reference.

Parameters:

z2d (np.ndarray) – 2D DEM array.
xres (float) – Cell size along x.
yres (float) – Cell size along y.
hillshade_az (float) – Hillshade azimuth in degrees.
hillshade_alt (float) – Hillshade altitude in degrees.
hillshade_z (float) – Hillshade vertical exaggeration factor.
search_radius_cells (int) – Ray walk length in cells.
n_directions (int) – Number of azimuth rays per cell.
svf_opacity (float) – SVF blend opacity in [0, 1].
po_opacity (float) – PO blend opacity in [0, 1].

Returns:

The VAT composite in [0, 1].

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.msrm_ref(z2d, xres, yres, fmin, fmax, x=2.0, edge_mode='reflect')

Multi-Scale Relief Model (Orengo & Petrie, 2018) verbatim port of MSRM.ipynb.

Constants:

\[ \begin{align}\begin{aligned}rr = (|xres| + |yres|) / 2\\base_i = max((fmin - rr) / (2 rr),\, 0)\\base_n = max((fmax - rr) / (2 rr),\, 0)\\i = \lfloor base_i^{1/x} \rfloor\\n = \lceil base_n^{1/x} \rceil\end{aligned}\end{align} \]

The MSRM is the mean of the relief models LP(round(k^x)) - LP(round((k+1)^x)) for k in [i, n-1].

Parameters:

z2d (np.ndarray) – 2D DEM array.
xres (float) – Cell size along x.
yres (float) – Cell size along y.
fmin (float) – Minimum scale in metres.
fmax (float) – Maximum scale in metres.
x (float) – Scale-progression exponent (default 2.0).
edge_mode (str) – Boundary handling mode for the box-mean lowpass.

Returns:

The MSRM map.

Return type:

np.ndarray

src.tests.raster_grid_postproc_ref.msrm_index_formula_ref(xres, yres, fmin, fmax, x=2.0)

Compute the MSRM scale indices (i, n) and the radii list.

Parameters:

xres (float) – Cell size along x.
yres (float) – Cell size along y.
fmin (float) – Minimum scale in metres.
fmax (float) – Maximum scale in metres.
x (float) – Scale exponent (default 2.0).

Returns:

Tuple (i, n, radii) where radii is the list of round(k^x) for k in [i, n].

Return type:

tuple