This is a typical case of a problem which is clear enough physically speaking, but mathematically messy. Where rigorous results are folkloristically employed to achieve some result which, actually, would need much more care in deriving it... But presumably, mathematical details would not change the physical picture. Here the difference between theoretical physics and mathematical physics evidently shows up.
Actually, that uniqueness theorem is not properly used in the example you mention. As it stands in your initial statement, the uniqueness property holds when $\Omega$ is an open and bounded subset of $\mathbb R^n$ and $\varphi$ is continuous in $\overline{\Omega}= \Omega \cup \partial \Omega$ and it is $C^2(\Omega)$ satisfying Poisson's equation $\Delta \varphi = \rho$ in $\Omega$ itself.
(The proof of uniqueness is a trivial consequence of a celebrated theorem on harmonic functions $\phi$ in $\Omega$, i.e., functions verifying $\Delta \phi =0$ in $\Omega$, which are continuous in $\overline{\Omega}$. That theorem establishes that, if $\Omega$ is open and it closure is compact, $\max_{\overline{\Omega}} |\phi|$ is reached in a point of $\partial \Omega$. Thinking of $\phi$ as the difference of two solutions of Poisson's equation filling the same boundary conditions, the uniqueness property easily arises.)
When $\Omega$ is not bounded, as in the mentioned example, where $\Omega = \{(x,y,z) \in {\mathbb R}^3\:|\: z>0\}$, one must add further requirements on the behaviour of $\varphi$ for $||(x,y, z)|| \to +\infty$, and there are a number of possibilities.
However the mentioned example suffers for another problem. In the above mentioned uniqueness result, $\rho$ is continuous because $\Delta \varphi$ is.
In the considered example instead $\rho$ is singular, properly speaking is a Dirac delta. There are several possibilities to deal with this problem.
The simplest one is replacing the point charge with a given spherically symmetric distribution - with total charge $q$ and confined in a bounded spherical small region - and continuously vanishing at the boundary of that region. In the rest of my answer I assume it.
Another possibility, technically more complicated is to remove the point occupied by the charge from $\Omega$. In this case the uniqueness result cannot be exploited as it stands because $\Omega$ acquires another part of boundary, where the potential diverges. Other approaches based on Green's identities instead of the maximum principle could be implemented in that case.
In this case in $\Omega$ the charge is $q$ (that is the $\rho$ distribution you mention) with distance $d$ from $\partial \Omega$ and $\varphi=0$ on $\partial \Omega$ because $\partial \Omega$ is a grounded conducting plane. The value of $\varphi$ on $\partial \Omega$ is constant, we are free to assume that it is $0$. The true boundary condition here is that $\varphi$ attains on $\partial \Omega$ the same value it attains for $||(x,y, z)|| \to +\infty$.
Let us pass to consider the situation where two charges stay at the reciprocal distance $2d$ along the $z$ axis, focussing on what happens in $\overline{\Omega}$ (not outside it) regarding charge distributions and boundary conditions of $\varphi$.
The $\rho$ distribution in the half space $\Omega = \{(x,y,z) \in {\mathbb R}^3\:|\: z>0\}$ is the same as in the previous case: There is the charge $q$ at distance $d$ from the plane at $z=0$.
Also the boundary conditions of $\varphi$ on $\partial \Omega$ and for $||(x,y, z)|| \to +\infty$ are the same as for the other case: The plane at $z=0$ is equipotential in view of the symmetry of the problem and the value of $\varphi$ thereon is the same as the value of $\varphi$ for $||(x,y, z)|| \to +\infty$.
Therefore, applying the uniqueness property in $\overline{\Omega}$, we are committed to conclude that the potential $\varphi$ in the region $\Omega \cup \partial \Omega$ is the same in both cases.
ADDENDUM. Actually one can use an argument arising from the theory of elliptic regularity to deal with uniqueness in situations where in a bounded region are present point charges described by Dirac deltas. The idea relies upon the following elliptic regularity result.
If $\phi$ is a distribution verifying $\Delta \phi =f$ (in weak sense) for a smooth ($C^\infty$) function $f$, then $\phi$ is a $C^\infty$ function up to zero measure set.
(It is worth stressing the the above result immediately entails the fantastic fact that harmonic functions are always $C^\infty$ and not only $C^2$, actually it is possible to prove that they are real analytic.)
This result leads to the following uniqueness theorem which can be improved making weaker some hypotheses on the behaviour of the function on the "regular" boundary.
Theorem. Suppose $\Omega \subset \mathbb R^n$ is non-empty open and $\overline{\Omega}$ is compact. Let $p \in \Omega$ and consider the problem:
$$\Delta \varphi(x) =\rho \quad x \in \Omega \setminus \{p\}$$
with boundary conditions
$$\varphi|_{\partial \Omega} = f$$
where
$$\varphi \in C^2(\Omega \setminus \{p\}) \cap C^0(\partial \Omega \cup \Omega \setminus \{p\} )$$ and $f \in C^0(\partial \Omega)$ and $\rho \in C^0(\Omega \setminus \{p\})$ are assigned.
If both $\varphi_1$ and $\varphi_2$ are solutions of the problem and $$ \lim_{x\to p} (\varphi_1(x)- \varphi_2(x))=0\:,$$
(where the limits can diverge or not exist, if considered separately in order to embody the case of a point charge at $q$) then
$$\varphi_1 = \varphi_2\:.$$
PROOF. With the given hypotheses, evidently $\phi:= \varphi_1-\varphi_2$ is continuous on $\Omega$, therefore it is a distribution for test functions, $h\in C_0^\infty (\Omega)$. If $B_\epsilon$ is a small ball around $p$ with radius $\epsilon$, using continuity of $\phi$ in particular, integrating by parts and defining $\Omega_\epsilon := \Omega \setminus B_\epsilon$, we have
$$\int_\Omega \phi \Delta h d^nx = \lim_{\epsilon \to 0^+}\int_{\Omega_\epsilon} \phi \Delta h d^nx =
\lim_{\epsilon \to 0^+} \int_{\Omega_\epsilon} (\Delta \phi) h d^nx =
\lim_{\epsilon \to 0^+} \int_{\Omega_\epsilon}(\rho-\rho) h d^nx$$ $$= \lim_{\epsilon \to 0^+} 0 = 0\:.$$
All that means that $\phi$ is a distribution solving $\Delta \phi =0$ in distributional sense. Therefore, in view of the mentioned elliptic regularity property, it is a smooth function up to a zero measure set. Since $\phi$ is continuous on $\Omega \setminus \{p\}$ and extends to a continuous function at $p$ (which has zero measure), $\phi = \varphi_1-\varphi_2$ is a smooth function everywhere in $\Omega$. In particular, by continuity of second derivatives the smoothly extended function $\phi$ verifies $\Delta \phi =0$ in the whole set $\Omega$ in the proper sense. By construction, we end up with a function $\phi$ which is $C^\infty(\Omega) \cup C^0(\overline{\Omega})$, satisfying $\Delta \phi =0$ in $\Omega$ and $\phi =0$ on $\partial \Omega$. In view of the standard uniqueness result, $\phi=0$ in $\overline{\Omega}$, i.e. $\varphi_1= \varphi_2$ in $\overline{\Omega}$. QED
