Spherical harmonics are not themselves full atomic orbitals. Consider the Hydrogen wave function, which separates into a radial part and an angular part. The latter is a spherical harmonic, but the former is some other function (in the case of Hydrogen it's a Laguerre polynomial). In general, we can approximate the angular part for other atoms with the same spherical harmonics, but we usually don't know the radial part analytically. Hence, if you want to directly evaluate Eq. (1) by integration you also need to find the radial part somewhere. The strength of Slater's and Koster's work, however, is that we can avoid this problem entirely, by hiding all radial dependence in the Slater-Koster parameters (in your example they are $V_{pd\sigma}$ and $V_{pd\pi}$).
For simplicity, I'll focus on the case of an $s$-$p$ overlap. This simplifies the geometry, and will let me borrow pictures and notation from this thesis. Let's say we have an $s$ orbital at site $i$ with wave function $\psi_{is}$, and a $p_\alpha$ orbital at site $j$ with wave function $\psi_{jp_\alpha}$, where $\alpha \in \{x,y,z\}$. In Dirac's bra-ket notation, the overlap between them can be written
$$
E_{s,p_\alpha}=\langle \psi_{is}|H_{2c}|\psi_{jp_\alpha}\rangle = \langle S|H_{2c}|P_\alpha\rangle, \tag{3}
$$
where $H_{2c}$ is the two-center Hamiltonian, and we've introduced short-hand notation for the the two wave functions. Of course, behind the bra-ket notation you have exactly the kind of overlap integral shown in Eq. (1).
Next step is to work out the geometry. Let $\vec{r}$ be the vector connecting sites $i$ and $j$, and let $\vec{d}$ be a unit vector along the same direction. We decompose the $p$ orbital at site $j$ into components parallel to ($\sigma$) and perpendicular to ($\pi$) the vector $\vec{d}$, as shown in this figure:
Figure from Anthony Carlson's 2006 MSc thesis at University of Minnesota.
In this notation, $\sigma$ and $\pi$ (and also $\delta$) is used to denote the component of angular momentum about the axis $\vec{d}$. $\sigma$ means zero, $\pi$ means $1$, etc. To proceed, we also define $\vec{a}$ as pointing along the $p$ orbital, and $\vec{n}$ as a vector perpendicular to $\vec{d}$ in the plane spanned by $\vec{a}$ and $\vec{d}$. Then, the $p$ orbital (at site $j$) can be decomposed
$$
|P_\alpha\rangle = \vec{a}\cdot\vec{d}|P_\sigma\rangle + \vec{a}\cdot\vec{n}|P_\pi\rangle. \tag{4}
$$
Then, the overlap is simply
$$
\langle S|H_{2c}|P_\alpha\rangle = \left( \vec{a}\cdot \vec{d} \right) \langle S | H_{2c} |P_\sigma\rangle + \left( \vec{a}\cdot\vec{n} \right) \langle S | H_{2c} | P_\pi\rangle, \tag{5}
$$
where the second term is zero by symmetry. Further, we can easily express this in terms of directional cosines since we can without loss of generality choose $\vec{a}$ to be parallel to one of the coordinate axes. With
$$
d_x=\frac{\vec{r}\cdot\hat{x}}{|\vec{r}|},\quad d_y=\frac{\vec{r}\cdot\hat{y}}{|\vec{r}|},\quad d_z=\frac{\vec{r}\cdot\hat{z}}{|\vec{r}|}, \tag{6}
$$
we get
$$
\langle S|H_{2c}|P_\alpha\rangle = d_\alpha \langle S|H_{2c}|P_\alpha\rangle = d_\alpha V_{sp\sigma}, \tag{7}
$$
where all radial overlap is hidden in the parameter $V_{sp\sigma}$. In tight-binding models it's common to denote this overlap integral $t_{sp\sigma}$ if it occurs in a hopping term.
The same approach works for other orbital combinations. You just need to set up the geometry and coordinate systems properly, and know where the atomic orbitals point. (Admittedly this can become quite complicated in some systems, e.g. transition metal oxides.) Then the Slater-Koster parameters can be treated as tuning parameters - either tuned to explore possible phenomena in some system, or fit to reproduce some experiment or calculated band structure.
