I can hopefully help you to see the fault in your logic by using some "physicist's thinking" — not very rigorous, but perhaps more illustrative (and hopefully the real mathematicians here won't tear me apart for this).
Here's a super-rough sketch by hand:
This is a sketch of one of your "slices" enclosed between two very close planes in heights $z$ and $z + \mathrm{d}z$, where by $\mathrm{d}z$ I mean an infinitesimally small increment of $z$. The distance of the curve from the axis of rotation is $r$ in the lower plane and slightly different $r + \mathrm{d} r$ in the upper one.
So what happens when you calculate the volume of the slice? You just calculate the volume of a small cylinder with radius $r$ and height $\mathrm{d} z$, getting $\pi r^2\,\mathrm{d} z$. Notice that in the sketch, this is the rectangle with two black horizontal and two blue vertical lines. So by counting like that, we're essentially "throwing away" the little colorful triangles.
Why do we get the correct answer? It's because the triangles are way too small and they give zero volume when rotated. You can intuitively see that just by looking at the areas in the sketch: you count the volume that corresponds to the black/blue rectangle with area $r\,\mathrm{d} z$, and throw away the triangle with area $\tfrac12 \mathrm{d} r\,\mathrm{d}z$. So the volume of the slice is infinitesimally small, and the volume of the triangle is infinitesimally small with respect to that infinitesimally small slice. It's "infinitesimally small squared"; no wonder it doesn't contribute.
However: when you calculate the surface area with your method, you're still counting the surface area of the little cylinders (i. e. the blue bits in the sketch). Sadly, that won't do, because the "real surface" is the orange one. And you're using $\mathrm{d} z$ where you should really use the length of the sloped orange line, which is $\sqrt{\mathrm{d} z^2 + \mathrm{d} r^2}$ (just apply the Pythagorean theorem to the colorful triangle). $\mathrm{d} r$ and $\mathrm{d} z$ are pretty similar in magnitude, so you can't just throw one of them out and expect to get away with it — it's like saying that $\sqrt2$ is pretty much equal to $1$ because $\sqrt2 = \sqrt{1 + 1}$ and we can neglect one of the $1$'s!
So, taking your half-sphere $r = \sqrt{R^2 - z^2}$, we first need to calculate
$$ \mathrm{d} r = \frac{\mathrm d (R^2 - z^2)}{2\sqrt{R^2-z^2}} = \frac{- z\,\mathrm d z}{\sqrt{R^2 - z^2}}. $$
And the surface area of the slice is essentially the surface of a cylinder with radius $r$ and height $\sqrt{\mathrm d z^2 + \mathrm d r^2}$, so it needs to be
$$ \mathrm{d} S = 2\pi r \sqrt{\mathrm d z^2 + \mathrm d r^2} = 2\pi \sqrt{R^2 - z^2} \times \sqrt{\mathrm d z^2 + \frac{z^2}{R^2 - z^2} \mathrm d z^2} = 2\pi R |\mathrm d z|. $$
Integrate this from $z = 0$ to $z = R$ to get the $2\pi R^2$ that we know is correct.
Best Answer
Mass
Imagine a lamina with constant density, $\rho$. The mass $m$ of the lamina and the area $A$ occupied by the lamina are in equal proportion. $1/4$ of the lamina has $1/4$ of the mass of the lamina, and so on. Now, what happens when the density isn't constant? This is what calculus was built for.
Say the density changes from point to point within a certain lamina. We can say that at a position $(x,y)$, the density $\rho(x,y)$ can be expressed as
$$\rho(x,y)=\frac{dm}{dA}$$
Now, to find the total mass of the object, we can integrate density with respect to volume:
$$dm = \rho(x,y) \ dA$$
$$m= \iint_R \rho(x,y) \ dA$$
So now that we have mass covered, let us define moments.
Moments
Say I give you a bat and a kettlebell. You're going to hold the bat horizontally and I will hang the kettlebell on different parts of the bat. In other words, I'm going to change the mass distribution of the system (bat and kettlebell) several times. You will find that the closer I hang the kettlebell to where you're holding, the easier the bat is to move. Mathematically speaking, suppose your hands are located at the origin. When I hang the kettlebell closer to your hands, the center of mass of the system moves closer to the origin. When I hang it further away from your hands, the center of mass moves farther from the origin. We turn the concept of how easily the bat moves into the idea of a moment.
We know that the further the kettlebell is from your hands, the harder the bat is to move. This means that moments are proportional to distance from the origin. Moments are proportional to mass, too: if I give you a $10$ kg kettlebell versus a $20$ kg kettlebell, you will undoubtedly notice a difference.
In 2D space, we use the $x$ and $y$ axes, so we will have moments relative to the $x$ axis and the $y$ axis. In other words, we will be rotating the bat around the two axes. Given the empirical evidence we have discussed, we can define the moments about the $x$ and $y$ axes as
$$M_x = m_1 y_1+m_2 y_2+...+m_n y_n = \sum_{i=1}^{n} m_i y_i$$
$$M_y = m_1 x_1 + m_2 x_2 + ... + m_n x_n = \sum_{i=1}^{n} m_i x_i$$
Now, the problem becomes complicated when instead of macroscopic kettlebells, you have microscopic changes in mass from point to point. Our friend calculus will be helpful here. Each piece of mass will be tiny, so our moments will become
$$M_x = \lim_{n \to \infty} y_1 \Delta m + y_2 \Delta m+...+ y_n \Delta m = \lim_{n \to \infty} \sum_{i=1}^{n} y_i \Delta m$$
$$M_y = \lim_{n \to \infty} x_1 \Delta m + x_2 \Delta m+...+ x_n \Delta m = \lim_{n \to \infty} \sum_{i=1}^{n} x_i \Delta m$$
Though we know that $\Delta m$ is just the density, $\rho (x,y)$, times the tiny piece of area $\Delta A$, so we now have
$$M_x = \lim_{n \to \infty} \sum_{i=1}^{n} y_i \Delta m = \lim_{n \to \infty} \sum_{i=1}^{n} y_i \rho (x,y) \Delta A = \iint_R y \rho (x,y) \ dA$$
$$M_y = \lim_{n \to \infty} \sum_{i=1}^{n} x_i \Delta m = \lim_{n \to \infty} \sum_{i=1}^{n} x_i \rho (x,y) \Delta A = = \iint_R x \rho (x,y) \ dA$$
Center of Mass
To find center of mass on any axis, we simply divide the moment about the opposing axis by the total mass. This means that
$$\overline{x} = \frac{M_y}{m} = \frac{\iint_R x \rho (x,y) \ dA}{\iint_R \rho(x,y) \ dA}$$
$$\overline{y} = \frac{M_x}{m} = \frac{\iint_R y \rho (x,y) \ dA}{\iint_R \rho(x,y) \ dA}$$
Hopefully this helps aid your understanding. I'm not too good with Riemann sums so I may need to edit this answer, but I think this should give you what you need.