As explained in the comments, the issue is not due to the presence of the charge. By superposition with the usual virtual charge solution, you are back to the case with no charges. More precisely, the real problem is a Laplace problem with $V(z=0)=1$ and $V(\infty)=0$ (normalising $V_0=1$).
You can already anticipate the problem. Due to invariance along the $x,y$ directions, the solution should be invariant along these directions as well. In order to match the boundary conditions at infinity, $V(z>0)=0$.
You can start to build intuition by making a dimension finite, breaking translational invariance to escape from the previous inconsistency. For example, say that $x\in[-L/2,L/2]$ with $V(x=-L/2,L/2)=0$ as well. You should recover the original problem when $L\to\infty$. By separation of variables (and translation invariance along $y$), you get:
$$
V = \frac{4}{\pi}\sum_{n=0}^\infty\frac{\sin[(2n+1)\pi (x/L+1/2)]e^{-(2n+1)\pi z/L}}{2n+1}
$$
The limit case $L\to \infty$ is not obvious from this expression as the simple limit is not summable. In this case, the fastest way to get the correct asymptotic behaviour is to express the sum in close form, using:
$$
\sum_{n=0}^\infty \frac{Z^{2n+1}}{2n+1} =\frac{1}{2} \ln\left(\frac{1+Z}{1-Z}\right)
$$
You get with $Z = ie^{\pi(ix-z)/L}$:
$$
V = \frac{4}{\pi}\Im\left[\frac{1}{2} \ln\left(\frac{1+Z}{1-Z}\right)\right]
$$
The branch of the logarithm is chosen from the boundary condition, i.e. when $Z = ie^{i\pi x}$, $V(Z)=1$ (geometrically, along the semicircle linking the two singularities at $Z=\pm1$). In the limit $L\to \infty$:
$$Z\to i$$
You therefore get the trivial limit as expected, since it is equivalent to fixing $L$ and sending $x,z\to 0$:
$$
V\to 1
$$
The solution is trivially harmonic, however, the boundary conditions at infinity did not "survive" the limit process. Mathematically, this is merely due to the fact that the convergence is not uniform on the entire domain. There is another way to understand this limit. From scaling invariance, fixing $x,y,z$ and sending $L\to\infty$ is the same as fixing $L$ and sending $x,y,z\to 0$. You are therefore sent to the boundary, and cannot detect the grounded boundary conditions.
One way to get a non-trivial limit would be to rescale $z$ as well, $z=L\tilde z$. Intuitively, this amounts to looking far from the grounded plane. In this case,
$$
Z\to ie^{-\pi \tilde z}
$$
so you get:
$$
V\to \frac{4}{\pi}\arctan(e^{-\pi\tilde z})
$$
Note that the solution is not harmonic anymore due to the non-isotropic scaling and the same phenomenon is observed for the boundary conditions at $x\to\infty$. However, it captures the correct exponential decay in the $z$ direction.
You could do a similar analysis with an insulated, grounded hemisphere or cube (minus a face) and sending the characteristic length to infinity. The scaling argument still applies, and shows that sending $L\to\infty$ amounts to looking at these solutions with $x,y,z\to0$. You will only see the boundary condition at $z=0$, so $V\to1$. In particular, you will observe the similar "disappearing" of boundary conditions.
These kind of non uniform limits are expected and illustrate that there is no natural consistent way of setting the boundary condition to $0$ at infinity. By some careful rescaling, you could get some interesting limiting behaviours, which could be useful depending on the problem at hand.
Hope this helps.
