Making a shakemap-grade PGA grid from a single station feed

A few months back I was asked to stand up a local PGA display for a small array of MEMS accelerometers feeding a Raspberry Pi over serial. Three stations, 100 Hz, one event at M4.1 about 18 km away. The obvious move would have been to pull ShakeMap, but the whole point was that these were our own sensors and nobody had a map product for them yet. I wanted a contour of peak ground acceleration that updated within a minute of the event, not a research-grade product.

The trap I hit first: raw accelerometer counts are not acceleration, and the baseline drifts. If you detrend with a single mean subtraction you get a nice clean zero that is wrong by the end of a 60 s window, and integration to velocity walks off the page. What worked was two passes: a running median (window ~1 s, 100 samples) to kill the DC step, then a 4th-order Butterworth high-pass at 0.1 Hz. That is deliberately aggressive at the low end. I care about the 0.5–15 Hz band where these sensors are actually trustworthy; below that the noise floor of a cheap MEMS part is doing the talking.

from scipy.signal import butter, sosfiltfilt
import numpy as np

sos = butter(4, 0.1, btype='highpass', fs=100, output='sos')

def pga_g(counts, sens_mv_per_g, adc_counts_per_mv):
    a = counts / adc_counts_per_mv / sens_mv_per_g
    a = a - np.median(a[:100])
    a = sosfiltfilt(sos, a)
    return 1000.0 * np.max(np.abs(a))  # mg

Then the grid. Three stations is not a grid, so I fitted a simple log-distance attenuation: pgm = A * exp(-k*R) with A and k from the three residuals, and let the contour routine interpolate in log space. That is a strong assumption. It is also honest if I say so on the plot.

Where it fails, in order of how often it bit me: (1) filtfilt is non-causal, so it smears energy before the first arrival — fine for peak amplitude, useless for onset picking; (2) three stations will always overfit a two-parameter curve, and a fourth bad station flips the map; (3) clipping — one sensor hit full scale on a M4.1 at 18 km and the peak was a lie, so I added a hard flag when counts exceed 95% of range; (4) the running-median detrend leaves a small residual step at the median window boundary that shows up as a fake spike in the differentiated trace, which is why I use the high-pass and not differentiation.

I still don’t have a good answer for how to report uncertainty on a three-station interpolation without fooling people into reading contour lines as data. If you’ve built small-array maps, how do you draw the confidence band — dashes, a shaded buffer around the stations, or just refuse to contour beyond the convex hull?