Existence of a uniformly continuous function $g$ on $\mathbb{R}$ where $f = g$ a.e.?

Yes.

Let $\{\phi_n\}$ be a sequence of standard mollifiers, so $\phi_n$ is continuous, $\phi_n \ge 0$, and $\int \phi_n = 1$. I claim the sequence $\{f \ast \phi_n\}$ is uniformly equicontinuous. For let $\epsilon > 0$ and choose $\delta$ so small that $|h| < \delta$ implies $\|f - f_h\|_{L^\infty} < \epsilon$. Then for any such $h$ and any $x \in \mathbb{R}$ we have $$\begin{align*}|(f \ast \phi_n)(x) - (f \ast \phi_n)(x+h)| &= |((f - f_h) \ast \phi_n)(x)| \\ &\le \int |f(y) - f_h(y)| \phi_n(x-y)\,dy \\ &\le \|f - f_h\|_{L^\infty} \cdot \int \phi_n(x-y)\,dy \\ &< \epsilon \cdot 1.\end{align*}$$ Moreover, $\|f \ast \phi_n\|_\infty \le \|f\|_\infty$. Passing to a subsequence if necessary, we have $f \ast \phi_n \to f$ almost everywhere. By Arzela-Ascoli, on any compact interval $[a,b]$ there is a subsequence of $f \ast \phi_n$ converging uniformly to some continuous function $g$, which must equal $f$ almost everywhere on $[a,b]$. By taking larger and larger intervals, and noting that $g$ is well defined, we get $g$ defined on all of $\mathbb{R}$, still having $g$ continuous and $f=g$ almost everywhere.

Now since $g$ is continuous, we have $\sup_x |g(x) - g_h(x)| = \|f - f_h\|_{L^\infty} \to 0$ as $h \to 0$. So in fact $g$ is uniformly continuous.

This is easier when passed to some sort of weak formulation. By Lebesgue differentiation theorem, for almost every x, $\lim_{r\to0} \frac{1}{|B_r|} \int_{x+B_r} f=f(x)$. Replace each f(x) by the left hand side. We have modified f on a set of measure 0, so the left hand side doesn't change further, and we have this equality everywhere. Also, the assumption $\lim_{h\to0} ||f_h-f||_\infty=0$ remains valid.

Now we show f is uniformly continuous. For $\epsilon>0$, there is $\delta>0$ such that for all $h<\delta$, $||f_h-f||<\epsilon$. Let $x-x'<\delta/3$. There is $r, r'<\delta/3$ such that $|f(x)-\frac{1}{|B_r|} \int_{x+B_r} f|<\epsilon$ and $|f(x')-\frac{1}{|B_{r'}|} \int_{x+B_{r'}} f|<\epsilon'$. Then $|f(x)-f(x')|<2\epsilon+\frac{1}{|B_r||B_{r'}|} \int_{y\in x+B_r} \int_{z\in x+B_{r'}} |f(y)-f(z)|=2\epsilon+\frac{1}{|B_r||B_{r'}|}\int_{h<\delta} \int_{y\in x+B_r, y+h\in x'+B_{r'}} |f(y)-f(y+h)|<2\epsilon+\frac{1}{|B_r||B_{r'}|}\int_{h<\delta} \int_{y\in x+B_r, y+h\in x'+B_{r'}} \epsilon=2\epsilon+\frac{1}{|B_r||B_{r'}|} \int_{y\in x+B_r} \int_{z\in x'+B_{r'}} \epsilon=3\epsilon$.