Kriging (NNGP) interpolation layer (CONNECTIONTYPE kriging)

[hlherrera]

Summary

Adds a new interpolation layer type, CONNECTIONTYPE kriging, implementing ordinary kriging as a localized Nearest-Neighbor Gaussian Process (NNGP). Like the existing IDW and Kernel-Density layers, it turns a point/vector source layer into a Raster surface and it additionally produces a second band with the Kriging standard error (predictive uncertainty), which IDW cannot.

Kriging and Gaussian-Process regression are the same method; the NNGP formulation (Datta et al., 2016, after Vecchia, 1988) predicts each output pixel from only its m nearest samples, so per-pixel cost scales with the neighbour count rather than the total sample count, it stays fast as the sample set grows.

How it works

• Slots into the existing interpolation framework next to IDW / kernel density (src/interpolation.c), fed by a source vector layer through CONNECTION (= the source layer's NAME), with the sample value taken from that layer's STYLE SIZE [attribute] exactly as IDW does.
• Numeric core in src/kriging_nngp.h (pure C, libc only): a uniform-grid spatial index for exact k-NN, a small dense LU solve per neighbourhood, and the ordinary/simple kriging predictor + variance.
• Emits two GDT_Byte bands 1 = predictive mean, 2 = predictive std dev selected with PROCESSING "BANDS=n" (default band 1).

PROCESSING keys

key	meaning	default
KRIGING_MODEL	EXPONENTIAL (Matérn) · GAUSSIAN · SPHERICAL	EXPONENTIAL
KRIGING_TYPE	ORDINARY (unknown local mean) · SIMPLE	ORDINARY
KRIGING_NEIGHBORS	nearest samples per pixel (the NNGP m)	16
KRIGING_RANGE	practical range, in pixels	auto
KRIGING_SILL	partial sill (value variance)	auto (sample variance)
KRIGING_NUGGET	nugget; 0 = exact interpolation, >0 = smoothing	0

Commits

1. kernel + framework wiring kriging.c, kriging_nngp.h, the MS_KRIGING connectiontype, dispatch + an N-band output path in interpolation.c, the == MS_IDW guards mirrored in mapraster.c / maprasterquery.c, and CMake.
2. lexer token adds kriging to maplexer.l and regenerates maplexer.c with the project's invocation flex --nounistd -Pmsyy -i. The large maplexer.c diff is the usual flex DFA-table reshuffle from adding a keyword (this file is excluded from clang-format in .pre-commit-config.yaml).
3. msautotest msautotest/gdal/kriging.map over a small, self-contained 16-point CSV, exercising the connectiontype, both bands, and the PROCESSING knobs.

Testing

• Builds clean against GDAL 3.13 / GEOS / PROJ using the project's Alpine CI recipe; kriging.c and the wiring compile with no warnings, and the new sources pass clang-format 15.0.7.
• The NNGP core has a standalone unit test: exact k-NN vs brute force, exactness at data points (zero variance there), and unbiasedness on a constant field.
• msautotest/gdal/kriging.map renders and matches the committed expected images. The surface is drawn by the vendored AGG with no fonts, so the expected PNGs are platform-independent.

Notes for reviewers

• kriging is now a reserved mapfile keyword (like idw), so a layer literally named kriging would need quoting.
• A fuller worked example (NZ rainfall, IDW vs kriging side by side, incl. the uncertainty band) is kept out of this PR because it depends on a larger sample dataset; happy to add a trimmed-down example if that would help.

[geographika]

Thanks, @hlherrera. The more generic parts of this PR seem good to me. I've no idea for the src/kriging.c code though. If AI assistance was used here, do you fully understand the code in case there are bugs/fixes in the future?

To add this feature you'll also need to add a docs page to https://github.com/MapServer/MapServer-documentation/.
It would be great to include the fully worked example you mentioned (A fuller worked example (NZ rainfall, IDW vs kriging side by side, incl. the uncertainty band)). It can be similar in structure to https://mapserver.org/output/idw.html

[hlherrera]

Thanks @geographika.

AI helped, but I worked through the implementation myself and I'm comfortable maintaining it. It's small: a dependency-free numeric core (kriging_nngp.h) with a unit test, plus thin MapServer glue (kriging.c) for the per-pixel solve over the m nearest samples. The msautotest case compares against committed expected images for regression coverage.

Docs: opened MapServer-documentation#1093, same style as the IDW page, with the NZ rainfall example and uncertainty band.