Green's function of 2D Poisson equation in the lower half-plane with Neumann boundary The solution of the 2D Poisson equation:
$$
\frac{\partial^2 V}{\partial x^2}+\frac{\partial^2 V}{\partial y^2}=\delta(x-x_s)\delta(y-y_s) \tag{*}
$$
subjected to the boundary condition $V(\infty)=0$ is:
$$
G(x,y)=\frac{1}{2\pi}\ln|r-r_s|
$$
Question: What is the solution of (*) subjected to the following boundary conditions?
$$
V(x\to\pm\infty)=0\\
V(y\to-\infty)=0\\
\frac{\partial V}{\partial y}\left|_{y=0}\right.=0
$$
 A: Since there's still confusion, I'll just write down the solution here
$$ V = \frac{1}{2\pi}\left(\ln\sqrt{(x-x_s)^2+(y-y_s)^2} + \ln\sqrt{(x-x_s)^2+(y+y_s)^2}\right) $$
The first term is the Green's function in free space. The second term is the same Green's function but its point source is $(x_s,-y_s)$
This page gives a more thorough explanation of the method of images (and also how to apply the Green's function). The half-plane example shown there is for a Dirichlet boundary, but you can easily infer a Neumann boundary by flipping the sign.
Key points:


*

*The reflected Green's function $G^*=G(x-x_s,y+y_s)$ is harmonic in the domain of interest, since its point source is outside the domain.

*The difference of $G-G^*$ is odd across the boundary, and therefore is Dirichlet.

*The sum of $G+G^*$ is even across the boundary, and therefore is Neumann.
A: I started working to this answer after seeing the comment of Dylan and the following comment of A Slow Learner, therefeore I was a little misplaced by the answer by Dylan himself (which however focuses more on the image method for constructing Green's functions): however I think it could be useful to see how the answer is obtained, and so I've decided to provide a(n almost) complete development below.
The Green's function $\mathscr{G}_N(x-x_s,y-y_s)$ for Neumann's problem for the Poisson equation (*) in the lower half plane is the solution of the following Dirichlet boundary value problem:
$$
\begin{cases}
\Delta \mathscr{G}_N(x-x_s,y-y_s)=\delta(x-x_s,y-y_s),\quad (x,y),(x_s,y_s)\in \mathbb{R}\times\mathbb{R}^-\\
\frac{\partial\mathscr{G}_N}{\partial y}(x-x_s,-y_s)=0
\end{cases}.\tag{1}\label{1}
$$
Note that $\frac{\partial\mathscr{G}_N}{\partial\mathbf{n}}=\frac{\partial\mathscr{G}}{\partial y}$ on the boundary $\{y=0\}$ of the lower half plane.
According to Vladimirov ((1983) §29.1, p. 368)), the general solution of problem \eqref{1} has the following form
$$
\mathscr{G}_N(x-x_s,y-y_s)=\frac{1}{2\pi}\ln|(x,y)-(x_s,y_s)|+g(x-x_s,y-y_s)
$$
where the first term on the right side is the fundamental solution of the laplacian while the second one is a function harmonic in the whole $\mathbb{R}\times\mathbb{R}^-$. In the problem posed by the OP,
$$
\mathscr{G}_N(x-x_s,y-y_s)=\frac{1}{2\pi}\big(\ln|(x,y)-(x_s,y_s)|+\ln|(x,y)-(x_s,-y_s)|\big).\tag{2}\label{2}
$$
Since $\frac{\partial\mathscr{G}_N}{\partial\mathbf{n}}=\frac{\partial\mathscr{G}_N}{\partial y}$ for problem \eqref{1}, by nothing that
$$
\begin{split}
\frac{\partial}{\partial y}\frac{1}{2\pi}\ln|(x,y)-(x_s,y_s)|& =\frac{\partial}{\partial y}\frac{1}{2\pi}\ln\sqrt{(x-x_s)^2+(y-y_s)^2}\\
&=\frac{1}{\pi}\frac{y-y_s}{(x-x_s)^2+(y-y_s)^2},
\end{split}
$$
it is easy to verify that the function \eqref{2} solves problem \eqref{1} and , by Green's formula it solves the general Neumann problem for the 2D Poisson's equation on the half plane.
[1] Vladimirov, V. S.  (1983)[1970], Equations of mathematical physics, Moscow: Mir Publishers, 2nd ed., pp. 464, MR0764399, Zbl 0207.09101 
(the Zbmath review refers to the first English edition).
