Observe that if you replace $P$ by $\bar{P} = P + ax + by + c$ and $S$ by $\bar{S} = S + bx - ay + d$, we have that
$$ \partial_x \bar{P} + \partial_y \bar{S} = \partial_x P + a + \partial_y S - a = \partial_x P + \partial_y S $$
and
$$ \partial_y \bar{P} - \partial_x \bar{S} = \partial_y P + b - \partial_x S - b = \partial_y P - \partial_x S $$
so strictly speaking, your source was wrong about uniqueness of the decomposition.
In more generality, the decomposition you mentioned is a manifestation of the Hodge decomposition, and in the case of three dimensions instead of two, it often goes by the name of Helmholtz decomposition.
Roughly speaking, we can compute $P$ by solving
$$ \triangle P = \operatorname{div} v = \partial_x v_x + \partial_y v_y $$
and
$$ \triangle S = \operatorname{curl} v = -\partial_x v_y + \partial_y v_x $$
That the decomposition is true hinges on
- The existence of of $P,S$ in a suitable function space solving the system of equations above (good candidate functions)
- The condition that in a suitable function space $\operatorname{div} v = 0$ and $\operatorname{curl}v = 0$ together implies that $v = 0$. (So that the candidate functions $P,S$ actually are the correct potentials; this will also imply uniqueness.)
This is the case, for example, if we assume that the derivatives of the coefficients $v_x$ and $v_y$ decays rapidly (faster than any polynomial). In that case using the Newton potential we can solve the Poisson equations for $P$ and $S$, and verify that they are up to constants the unique bounded potentials for the decomposition (notice that $ax + by + c$ is not bounded unless $a = b = 0$).