I'm working with Griffiths Electrodynamics, and he introduces a uniqueness theorem:
First Uniqueness Theorem: The potential V in a volume Ω is uniquely determined if (a) the charge density throughout the region, and (b) the value of V on the boundary ∂Ω, are specified.
I'm a little confused how Griffiths uses this theorem in examples: In the classic image problem (finding the potential due to a point charge q a distance d above an infinite, grounded conducting plane), the trick is to forge the original problem and the configuration of q and its mirror image of the point charge q through the plane (and this new charge is −q). The boundary conditions are given (V=0 on the plane, and V→0 far from the point charge), but how do we know that the charge distribution ρ is the same in this second scenario as it is in the first?
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 Ω is an open and bounded subset of Rn and φ is continuous in ¯Ω=Ω∪∂Ω and it is C2(Ω) satisfying Poisson's equation Δφ=ρ in Ω itself.
(The proof of uniqueness is a trivial consequence of a celebrated theorem on harmonic functions ϕ in Ω, i.e., functions verifying Δϕ=0 in Ω, which are continuous in ¯Ω. That theorem establishes that, if Ω is open and it closure is compact, max 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