Yep, $f$ is convex. Indeed, even the assumption that $f \in C^1$ can be dispensed with, and replaced with continuity and a local subgradient condition. That is, if, for every $x \in \Bbb{R}^n$, there exists some $v_x \in \Bbb{R}^n$ and $\varepsilon_x > 0$ such that
$$\|y - x\| < \varepsilon_x \implies f(y) \ge f(x) + v_x^\top (y - x), \tag{$\star$}$$
then $f$ is convex.
We show that each local subgradient can be replaced by a global subgradient. That is, we may remove the condition $\|y - x\| < \varepsilon_x$ in $(\star)$. Fix $x_0 \in \Bbb{R}^n$ and consider the function
$$g(y) = f(y) - v_{x_0}^\top(y - x_0) - f(x_0).$$
Note that $g$ satisfies the condition $(\star)$ as well, $g(x_0) = 0$, and $g$ has a local minimum at $x_0$.
Fix $y_0 \in \Bbb{R}^n \setminus \{x_0\}$. We wish to show that the conclusion of $(\star)$ holds when $y = y_0$ and $x = x_0$, even if we don't have $\|y_0 - x_0\| < \varepsilon$. In particular, we are showing that $g(x_0)$ is a global minimum, not just a local one.
Restricting $g$ to the line segment $[x_0, y_0]$, we know that $g$ is continuous, and hence it attains a global maximum on this interval. This maximum must be at least $0$, since $g(x_0) = 0$. Further, since $x_0$ is a local minimum of $g$, this maximum must be achieved somewhere other than $x_0$.
Pick a point $x_1 \in \operatorname{argmax}_{x \in [x_0, y_0]} g(x) \setminus \{x_0\}$. If $y_0$ lies in this set, then we are done, so assume this is not the case, and hence $x_1 \neq y_0$. As per $(\star)$, there exists some $w_{x_1} \in \Bbb{R}^n$ and $\varepsilon_{x_1} > 0$ such that
$$\|y - x_1\| < \varepsilon_{x_1} \implies g(y) \ge g(x_1) + w_{x_1}^\top(y - x_1).$$
Note that, for sufficiently small $\lambda$, $z_\lambda := x_1 + \lambda(y_0 - x_0) \in [x_0, y_0]$ and $\|z_\lambda - x_1\| < \varepsilon$, hence
$$g(x_1) \ge g(z_\lambda) \ge g(x_1) + w_{x_1}^\top(z_\lambda - x_1),$$
which implies that
$$0 \ge w_{x_1}^\top(z_\lambda - x_1) = \lambda w_{x_1}^\top(y_0 - x_0)$$
for all sufficiently small $\lambda$. This implies that $w_{x_1}^\top(y_0 - x_0) = 0$, and hence, for sufficiently small $\lambda$,
$$g(x_1) \ge g(z_\lambda) \ge g(x_1) + 0,$$
hence $z_\lambda$ also maximises $g$. That is, we have $x_1$ is in the interior of $\operatorname{argmax}_{x \in [x_0, y_0]} g(x)$, relative to the line segment $[x_0, y_0]$, and hence this $\operatorname{argmax}$ set is open in $[x_0, y_0]$. On the other hand, given the continuity of $g$, it is also closed, and hence must be the full line segment. But, we explicitly assumed $y_0$ did not maximise $g$ on the line segment, hence we have a contradiction. Thus, the maximum of $g$ occurs at $y_0$, hence $g$ achieves a global minimum at $x_0$.
This means that $f$ has a subgradient $v_x^\top$ at every point $x$. This implies the convexity of $x$. In particular, we may express:
$$f(y) = \sup_{x \in \Bbb{R}^n} (f(x) + v_x^\top(y - x)),$$
which makes $f$ the pointwise supremum of affine functions, which is convex.