Not really, but it depends on what you consider as "explicit". For example, you could take the 2D Fourier transform of your problem and the solution will depend explicitly on the Fourier transforms of $f_1$ and $f_2$.
Nonetheless, I see a more direct trick. I don't know in which context your system of PDEs arises, but it can be seen as inhomogeneous Cauchy-Riemann equations for the complex variable $z = x+iy$ and the function $g(z,\bar{z}) = g_1(x,y) + ig_2(x,y)$, with the source term $f(z,\bar{z}) = \frac{1}{2}(f_1(x,y) + if_2(x,y))$, such that $\partial_{\bar{z}}g = f$, where $\partial_{\bar{z}} = \frac{1}{2}(\partial_x + i\partial_y)$ is the Wirtinger derivative. N.B. : $f$ and $g$ are not holomorphic in that case, that is why they also depend on the conjugate variable $\bar{z}$.
The solution is then given by
$$
g(z,\bar{z}) = A(z) + \int_D \frac{f(z',\bar{z}')}{z'-z} \frac{\mathrm{d}z'\wedge\mathrm{d}\bar{z}'}{2\pi i},
$$
where $D\subset\mathbb{R}^2$ is your said region and $A(z)$ an holomorphic function on $D$. So, the antiderivatives of $f_1,f_2$ do not appear in the solution, because it is basically made of the convolution of $f$ with the kernel $\frac{1}{2\pi iz}$.