$\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}}
\newcommand{\PP}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}
\newcommand{\tf}{\widetilde{f}}$
As you noted,
\begin{equation*}
\tf-f=\frac{1}{\ep}\sum_n a_n 1_{[n\ep,(n+1)\ep)},
\end{equation*}
where
\begin{equation*}
a_n(x):=\int_{n\ep}^{(n+1)\ep} \int_x^s f'(\tau) d \tau \ ds.
\end{equation*}
For $x\in[n\ep,(n+1)\ep)$,
\begin{equation*}
|a_n(x)|\le\int_{n\ep}^{(n+1)\ep}\int_{n\ep}^{(n+1)\ep} |f'(\tau)| d \tau\ ds
=\int_{n\ep}^{(n+1)\ep} |f'(\tau)| d \tau\ \ep\le b_n\ep\sqrt\ep,
\end{equation*}
where
\begin{equation*}
b_n:=\sqrt{\int_{n\ep}^{(n+1)\ep} |f'(\tau)|^2 d \tau}.
\end{equation*}
So,
\begin{equation*}
|\tf-f|\le\sum_n b_n \sqrt\ep \, 1_{[n\ep,(n+1)\ep)}
\end{equation*}
and
\begin{equation*}
\|\tf-f\|_2\le\sqrt{\sum_n b_n^2 \ep^2}=\ep\|f'\|_2\le\ep\|f\|_{W^{1,2}},
\end{equation*}
as desired.
Consider now the general case when $f$ is a function on $\R^d$ with $d\ge2$ and
\begin{equation*}
\tf=\frac{1}{\ep^d}\sum_{n\in\Z^d} 1_{\de_n}\int_{\de_n} ds\,f(s),
\end{equation*}
where $\de_n:=\prod_1^d\,(n_i\ep,(n_i+1)\ep)$ for $n=(n_1,\dots,n_d)$.
Take any real
\begin{equation*}
p>d
\end{equation*}
and let $q:=\frac p{p-1}$, so that $\frac1p+\frac1q=1$.
We have
\begin{equation*}
\tf-f=\frac{1}{\ep^d}\sum_n a_n 1_{\de_n}, \tag{!}
\end{equation*}
where
\begin{equation*}
a_n(x):=\int_{\de_n} ds \int_0^1 dt\, f'(x+t(s-x))\cdot(s-x)
\end{equation*}
$f':=\nabla f$ and $\cdot$ is the dot product.
By the change of variables from $s$ to $\tau=x+t(s-x)$, for $x\in\de_n$ we have
\begin{equation*}
a_n(x)=\int_0^1 dt\int_{\de_n} \frac{d\tau}{t^d} \, f'(\tau)\cdot\frac{\tau-x}t\,1_{(t_*,\infty)}(t),
\end{equation*}
where
\begin{align*}
t_*:=t_*(n,x,\tau)&:=\max_1^d\Big(\frac{\tau_i-x_i}{x_i-n_i\ep}\vee\frac{x_i-\tau_i}{(n_i+1)\ep-x_i}\Big) \\
&\ge\max_1^d\frac{|\tau_i-x_i|}{\ep}\gg\frac{|\tau-x|}{\ep},
\end{align*}
with the constants in $\gg$ and $\ll$ depending only on $d$ and $p$, and
$|\ \, |$ also denoting the Euclidean norm on $\R^d$. Note that
\begin{equation*}
\int_{t_*}^\infty\frac{dt}{t^{d+1}}\ll\frac1{t_*^d}\ll\frac{\ep^d}{|\tau-x|^d}.
\end{equation*}
It follows that for $x\in\de_n$
\begin{equation*}
|a_n(x)|\ll\ep^d\,\int_{\de_n} \frac{d\tau}{|\tau-x|^{d-1}} \, |f'(\tau)|
\le\ep^d\,b_n c^{1/q}, \tag{!!}
\end{equation*}
by H\"older's inequality, where
\begin{equation*}
b_n:=\Big(\int_{\de_n}d\tau\,|f'(\tau)|^p\Big)^{1/p},
\end{equation*}
\begin{equation*}
c:=\int_{|\tau-x|\le\ep\sqrt d}\frac{d\tau}{|\tau-x|^{(d-1)q}}
\ll\int_0^{\ep\sqrt d}\frac{dr}{r^{(d-1)q}}\,r^{d-1}
\ll\ep^{(p-d)/(p-1)}.
\end{equation*}
Thus, by (!) and (!!),
\begin{equation*}
\|\tf-f\|_p\ll\Big(\sum_n b_n^p \ep^{p-d}\ep^d\Big)^{1/p}=\ep\|f'\|_p\le\ep\|f\|_{W^{1,p}},
\end{equation*}
as desired.
