The surface is defined by $$z^2=-(x^2+y^2-16)((x-2)^2+y^2-1)((x+2)^2+y^2-1).$$
How would I get the total Gaussian curvature?
I'm aware of that if the surface is represented by $X=(x,y,z)$, then the Gaussian curvature is $$k=\frac{LM-N^2}{EG-F^2}$$ where $E=X_x\cdot X_x,\, F=X_x\cdot X_y,\, G=X_y\cdot X_y$ and $L=X_{xx}\cdot n,\, M= X_{xy}\cdot n,\, N=X_{yy}\cdot n$.
Here $$n=\frac{X_x\times X_y}{|X_x\times X_y|},$$ the unit normal vector.
But the calculation is extremely overwhelming, because of the square power of $z$. Is there an other way other than just calculating everything?
And my definition of the total curvature is the surface integral of Gaussian curvature $\int_S k\,dA$, where $S$ is the surface defined above. Is my definition correct?
Not clear to me what you want.
In case you want $\int KdA$
Your definition is OK, it implies evaluation for the entire surface. This is a topological constant or invariant, a part of Gauss Bonnet theorem aka Integral Curvature. It is one of constituents in the theorem connecting isometric invariants and topological invariants introduced in such a beautifully permanent way. Briefly
$$ \int k_g ds + \sum_{i=1}^{n}\psi_{i }+\int KdA= 2 \pi(1-\chi) $$
The first two are isometric invariants that vanish for continuous geodesics.
The topological invariant $ KdA $ also equals $ 2 \pi (2-2 g) $
where $\chi $ is the Euler characteristic. No need to make any calculations. The geometry /topology should be understood.
For a sphere $g=1$ and the total solid angle is $4 \pi$ steradians. Also $ \chi= V+F-E,$ for polyhedral surfaces, for details please google.
We have to look to the genus or how many non-intersecting holes the solid body has. The given surface is triply connected with two holes i.e., topologically it is a sphere with two holes $(g=2) $ drilled in. Mid-section shows the two holes.
There are three surfaces, a sphere and two eccentric cylinders.
$$ (x^2+y^2-16),\; (x-2)^2+y^2-1).\;(x+2)^2+y^2-1)$$
that are mathematically blended by multiplication in the equation you have given. So this a good example of a blended surface where two hollow tubes are blended into a sphere with a small degree of homeomorphic deformation.
Here there are two holes so integral curvature is
$$ 2 \pi \chi = 2 \pi (2-2g)= - 4\pi$$
In case you want Gauss curvature $K$ only
There is a formula suitable surfaces given in the Mongé form:
$$ z=f(x,y);\;$$
$(p,q)$ are first partial derivatives;$(r,t,s),$ second partial derivatives:
$$K=\dfrac{(rt-s^2)}{(1+p^2+q^2)^2}$$
Partially differentiate twice and plug in. Shall copy the humongous output by
Mathematicaif you want.