# Smooth RF coverage maps in WebGL: a production playbook

**Replace one line of shader code and your grid squares vanish.** The single highest-impact fix for pixelated RF coverage overlays is swapping hardware bilinear texture sampling for a **Catmull-Rom bicubic fragment shader**—9 texture fetches instead of 1, near-zero performance cost, and mathematically correct C1-continuous interpolation that passes through your actual RSRP values. Professional RF tools like Atoll and CloudRF sidestep the problem entirely by computing at every grid cell, but when you're rendering a sparse client-side grid (1,975–6,675 points), shader-based interpolation is the correct approach. This report covers exactly how to implement it, what the industry does, and which open-source libraries can accelerate the work.

## How professional RF tools avoid the problem you're solving

Professional RF planning tools don't interpolate sparse data—they brute-force compute signal values at every pixel. **CloudRF** runs its SLEIPNIR propagation engine server-side at user-specified resolution (down to 1 m with LiDAR), outputs GeoTIFF or pre-colored PNG, and the client simply overlays the image via `L.imageOverlay` in Leaflet or drapes it on CesiumJS for 3D. There is no client-side interpolation. **Atoll** (Forsk) computes predictions at **5–50 m grid resolution** within its desktop GIS—smoothness comes from grid density, not post-processing. **SPLAT!** does the same with Longley-Rice/ITM, outputting PPM rasters where each pixel maps 1:1 to a DEM grid cell; its `-sc` flag adds smooth color gradients but no spatial interpolation.

The pattern is universal: propagation model → dense raster → colorize → overlay. Crowdsourced coverage services like **Ookla/Speedtest** (which uses Mapbox GL JS with WebGL rendering) and **OpenSignal** follow a different path: they aggregate sparse measurements into spatial bins, apply kernel density estimation or IDW smoothing server-side, then serve the result as raster tiles. The industry-standard interchange format is **Cloud-Optimized GeoTIFF (COG)** with bilinear or cubic resampling via GDAL at tile-generation time.

Your situation is fundamentally different from both. You have a moderate-density regular grid arriving from the backend, and you need the client to render it smoothly at arbitrary zoom. This makes shader-based interpolation the right tool—not denser computation, not server-side tiling.

## Catmull-Rom in a fragment shader: the production solution

The core insight is that your current pipeline—upload grid as float texture, sample with `GL_LINEAR`, apply colormap—is already 90% correct. Hardware bilinear interpolation produces C0 continuity (values match at grid edges, but derivatives don't), causing visible seams. **Catmull-Rom spline interpolation** provides C1 continuity (smooth first derivatives) while still passing through your exact data values, unlike B-spline bicubic which smooths/blurs peaks.

The 9-tap Catmull-Rom implementation, widely used in production (ported from TheRealMJP's gist with 108 GitHub stars), replaces your `texture2D()` call:

```glsl
vec4 SampleTextureCatmullRom(sampler2D tex, vec2 uv, vec2 texSize) {
    vec2 samplePos = uv * texSize;
    vec2 texPos1 = floor(samplePos - 0.5) + 0.5;
    vec2 f = samplePos - texPos1;

    vec2 w0 = f * (-0.5 + f * (1.0 - 0.5 * f));
    vec2 w1 = 1.0 + f * f * (-2.5 + 1.5 * f);
    vec2 w2 = f * (0.5 + f * (2.0 - 1.5 * f));
    vec2 w3 = f * f * (-0.5 + 0.5 * f);

    vec2 w12 = w1 + w2;
    vec2 offset12 = w2 / (w1 + w2);

    vec2 texPos0  = (texPos1 - 1.0) / texSize;
    vec2 texPos3  = (texPos1 + 2.0) / texSize;
    vec2 texPos12 = (texPos1 + offset12) / texSize;

    vec4 result = vec4(0.0);
    result += texture2D(tex, vec2(texPos0.x,  texPos0.y))  * w0.x  * w0.y;
    result += texture2D(tex, vec2(texPos12.x, texPos0.y))  * w12.x * w0.y;
    result += texture2D(tex, vec2(texPos3.x,  texPos0.y))  * w3.x  * w0.y;
    result += texture2D(tex, vec2(texPos0.x,  texPos12.y)) * w0.x  * w12.y;
    result += texture2D(tex, vec2(texPos12.x, texPos12.y)) * w12.x * w12.y;
    result += texture2D(tex, vec2(texPos3.x,  texPos12.y)) * w3.x  * w12.y;
    result += texture2D(tex, vec2(texPos0.x,  texPos3.y))  * w0.x  * w3.y;
    result += texture2D(tex, vec2(texPos12.x, texPos3.y))  * w12.x * w3.y;
    result += texture2D(tex, vec2(texPos3.x,  texPos3.y))  * w3.x  * w3.y;
    return result;
}
```

**Performance is essentially free.** Benchmarks on a GTX 980 show 9-tap Catmull-Rom at 1920×1080 costs **~0.32 ms** versus ~0.30 ms for single-tap bilinear. Your ~80×85 texel coverage texture is trivial. The critical rule: **interpolate raw scalar RSRP values first, then apply the colormap**. If you interpolate after colorization, you get color-space artifacts (muddy intermediate colors between discrete bands).

For the absolute quickest improvement with zero extra texture fetches, Inigo Quilez's **smoothstep coordinate remapping trick** eliminates grid edges with a single line change:

```glsl
vec4 textureSmooth(sampler2D tex, vec2 uv, vec2 texSize) {
    vec2 p = uv * texSize + 0.5;
    vec2 i = floor(p);
    vec2 f = p - i;
    f = f * f * f * (f * (f * 6.0 - 15.0) + 10.0); // quintic hermite
    return texture2D(tex, (i + f - 0.5) / texSize);
}
```

This gives C2 continuity with a single texture read, though it introduces slight positional bias. For a quick visual test before implementing full Catmull-Rom, it's ideal.

## When to use IDW instead, and how the GPU handles it

Catmull-Rom works because your data sits on a **regular grid**. If your backend ever returns irregular point distributions (e.g., drive-test measurements, crowdsourced data), **Inverse Distance Weighting (IDW)** on the GPU becomes the right approach. The production technique uses WebGL's additive blending to parallelize across data points rather than looping in a shader:

For each data point, render a full-screen quad where the fragment shader computes `w = 1/dist^p` and outputs `(w × value, w, 0, 1)` into a framebuffer with `gl.blendFunc(gl.ONE, gl.ONE)`. After N passes (one per point), a final shader reads the accumulated texture and computes `interpolated_value = R_channel / G_channel`. This avoids the O(N) per-pixel loop that would choke a fragment shader at 7,000 points.

The open-source **`mapbox-gl-interpolate-heatmap`** library implements exactly this pattern, with a `framebufferFactor` parameter (typically 0.3–0.5) that renders the IDW computation at reduced resolution for performance, then upscales for display. The **`temperature-map-gl`** library from ham-systems provides the same approach in a minimal ~3 KB package.

**For your regular-grid case, IDW is strictly worse than Catmull-Rom**: it's slower (N draw calls vs. 9 texture fetches), creates bull's-eye artifacts around data points, and doesn't exploit grid structure. Reserve it for irregular data.

## Mapping library internals and the tiled raster alternative

If you want to move beyond a custom WebGL overlay, understanding how major mapping libraries handle this problem reveals useful architectural patterns:

**Mapbox GL JS** exposes `raster-resampling: 'linear'` (bilinear) or `'nearest'` for raster tile layers—standard GPU texture filtering, same limitation you're hitting. Its built-in heatmap layer uses Gaussian kernel density estimation with additive blending into an offscreen half-float texture, designed for point density visualization, **not scalar interpolation**. However, Mapbox v3's `raster-array` source with `raster-color` expressions enables client-side colorization of raw data tiles, which is directly applicable. **deck.gl** offers `BitmapLayer` with configurable `textureParameters` (set `magFilter: 'linear'`) and supports full custom shader injection via `fs:DECKGL_FILTER_COLOR` hooks, but its `HeatmapLayer` is again KDE-based, not interpolation-based.

For a **tiled architecture** (useful if your coverage areas grow large), the proven pipeline is:

- Backend computes RSRP grid → stores as GeoTIFF
- GDAL resamples with `cubicspline` or `lanczos`: `gdalwarp -r cubicspline -tr 10 10 input.tif upsampled.tif`
- `gdal2tiles.py` generates XYZ tile pyramid
- Client displays via `L.tileLayer` with standard bilinear filtering

The **IHME `leaflet.tilelayer.glcolorscale`** library offers a sophisticated variant: encode 32-bit float values into PNG RGBA channels server-side, decode in a WebGL fragment shader client-side, and apply dynamic color scales without re-tiling. Its companion `leaflet.tilelayer.gloperations` adds GPU-based convolution smoothing. This pattern preserves raw values for pixel queries while enabling dynamic color ramp changes.

## Handling boundaries, gaps, and color mapping

Three practical concerns beyond core interpolation deserve specific solutions. For **coverage boundaries**, apply `smoothstep` fading based on signal strength or a validity mask:

```glsl
float signalStrength = SampleTextureCatmullRom(u_data, uv, texSize).r;
float boundaryAlpha = smoothstep(-115.0, -105.0, signalStrength); // fade near noise floor
gl_FragColor = vec4(colormap(t), boundaryAlpha * u_opacity);
```

For **missing grid cells**, store a validity mask as a second texture channel (or use a sentinel value like -9999). Interpolate both the value and mask with Catmull-Rom; the interpolated mask naturally creates smooth alpha transitions at data boundaries without hard edges.

For **color mapping**, use a **1D texture lookup** rather than branching `if/else` chains in GLSL. Upload your RSRP→color ramp as a 256×1 RGBA texture, normalize your interpolated value to [0,1], and sample: `vec3 color = texture2D(u_colormap, vec2(t, 0.5)).rgb`. This is faster than arithmetic color functions and trivially supports any color ramp. The `glsl-colormap` package provides standard scientific palettes (viridis, jet, etc.) as pure GLSL functions if you prefer avoiding the extra texture.

## Recommended implementation path

The optimal architecture for your Electron + React + Leaflet stack, given 1,975–6,675 grid points:

**Immediate fix (30 minutes):** Replace `texture2D(u_data, uv)` with the smoothstep trick in your existing fragment shader. This eliminates visible grid squares with zero performance cost and zero architectural changes.

**Production implementation (1–2 days):** Implement the 9-tap Catmull-Rom `SampleTextureCatmullRom()` function. Pack your RSRP grid into a float texture (R channel = RSRP value, G channel = validity mask). Apply Catmull-Rom to both channels. Use a 1D colormap texture for color mapping. Add smoothstep boundary fading. This produces results visually indistinguishable from CloudRF-quality coverage maps.

**If you outgrow the single-texture approach** (very large coverage areas, multiple overlapping cells): transition to the float-encoded tile pipeline using `leaflet.tilelayer.glcolorscale` or generate server-side tiles with GDAL cubic resampling. The tile pyramid handles LOD automatically and scales to arbitrarily large coverage areas.

## Conclusion

The gap between your current output and professional RF coverage visualization isn't architectural—it's a single shader function. Professional tools achieve smoothness through brute-force grid density (computing at every cell), but shader-based Catmull-Rom interpolation produces equivalent visual quality from sparse grids at negligible GPU cost. The 9-tap implementation requires **9 texture fetches** versus your current 1, adds C1 continuity that eliminates all visible grid edges, and preserves exact RSRP values at grid points—unlike Gaussian blur or B-spline smoothing, which distort the data. For truly irregular point data, GPU-accelerated IDW via additive blending is proven in production libraries like `mapbox-gl-interpolate-heatmap`. The critical implementation principle: always interpolate raw scalar values first, colorize second.