Everything you say about centers of the linearized differential system vs centers of the true differential system holds. To show that a "linearized" center is a "true" center, I am afraid one is led to find an invariant of the dynamics, at least in a neighborhood of the center--and for that, I know no general approach.
In principle, some quantity $H(x,y)$ is invariant by the differential system
$$\dot x=f(x,y),\qquad\dot y=g(x,y),$$ if and only if $H$ solves the PDE
$$
\partial_xH(x,y)\cdot f(x,y)+\partial_yH(x,y)\cdot g(x,y)=0.
$$
hence to find an invariant of the dynamics is equivalent to solving this PDE.
In the present case however, finding an invariant of the dynamics is relatively direct, once one notes that some quantity
$$
H(x,y)=F(x)+G(y)
$$
is invariant by the differential system
$$
\dot x=f(x,y),\qquad\dot y=g(x,y),
$$
if and only if
$$
F'(x)\cdot f(x,y)+G'(y)\cdot g(x,y)=0.
$$
Here, one gets
$$
\frac{x}{d-cx}F'(x)=\frac{y}{a-by}G'(y),
$$
hence both sides are constant, say without loss of generality that they are both equal to $1$. Thus,
$$
F'(x)=\frac{d}x-c,\qquad G'(y)=\frac{a}y-b,
$$
which shows that an invariant of the dynamics on $xy\ne0$ is
$$
H(x,y)=d\log|x|+a\log|y|-cx-by.
$$
To get an invariant of the dynamics on the whole $(x,y)$-plane, consider
$$
H(x,y)=|x|^d\,|y|^a\,\mathrm e^{-cx-by}.
$$