The question seems to have been edited... I saw the query about $\sqrt{z}+(z-2)^{1/3}$. I'll do the revised version afterward. The branch points are where either of the two radicals behaves strangely (has fewer than the expected number of distinct roots), namely, 0 and 2.
Over every point other than 0, there are two branches of the square root. To make two single-valued square roots, "cutting" the plane along any ray from 0 to infinity works.
Over everypoint other than 2, there are three branches of that cube root. To make three single-valued cube roots, cut the plane along any ray from 2 to infinity.
To take care of both radicals simultaneously, basically we need two rays, one from 0 to infinity, the other from 2 to infinity. We could take a ray that starts at 0, passes through 2, and goes to infinity, but this is a slightly weird choice, and doesn't really help anything much.
The Riemann surface has to be/is a "ramified covering" of the plane on which the function is well-defined. This requires a six-sheeted covering to incorporate the $2\cdot 3$ sheets needed by the square root and the cube root.
To determine the genus via the Riemann-Hurwitz formula, observe that the point at infinity has six-fold ramification, giving one point with ramification degree 6. There are 3 different points with ramification degree 2, and 2 points with ramification degree 3. Thus, by Riemann-Hurwitz,
$$ 2 - 2g' = 6(2 - 2g) - \sum_P (e_P-1)
= 12 - ((6-1)+3(2-1)+2(3-1)) = 12 - 12$$
Thus, the genus of the cover is $g'=1$. [Edited: garbled Riemann-Hurwitz!]
With $\sqrt{z}+(z-4)^{1/4}$, there are 8 sheets altogether. The point at infinity behaves a little differently, since the cover on which $(z-4)^{1/4}$ exists near infinity already has a $\sqrt{z}$ there. Thus, there are two points lying over the point at infinity, each with ramification degree 4. There are 4 points lying over 0, each with ramification 2, and 2 points over 4, each of ramification 4. Riemann-Hurwitz is
$$
2-2g' = 8(2 - 2g) - \sum_P(e_P-1) = 16 - (2(4-1)+4(2-1)+2(4-1))
= 16 - 16$$
so the genus of the cover is $g'=1$.
The "branch-cuts" along the imaginary axis seem to be a software artifact. The value of $\sqrt{z^2-1}$ is (apparently) computed by first computing $z^2-1$, and then taking the square root, choosing the square root with non-negative real part [and probably the one with non-negative imaginary part when the square root is purely imaginary]. That loses the information where the value whose square root shall be taken comes from, and thus leads to an even function (since $z^2-1$ is even), instead of a continuous branch, which would be an odd function (with the branch-cut $[-1,1]$). To avoid that, you need to remember where $z$ lies, and for $z$ in the left [or in the right, if you want the other branch] half-plane multiply the computed value with $-1$ [also for $z$ in the upper or lower half of the imaginary axis].
Then you should get a nice smooth plot, with a branch-cut on the interval $[-1,1]$.
For branch cuts I thought we had to identify where the argument of the square root is zero, which gives us the branch points and then we are free to take lines from those branch points to infinity in any way we like.
That describes the situation for $\sqrt{z}$ and similar, but generally, there are a few other considerations. If you have a rational function $f$, and consider $\sqrt{f(z)}$, the potential branch points are the zeros and poles of $f$.
But zeros and poles of even multiplicity are not actual branch points (in a neighbourhood of such a point you have $f(z) = (z-z_0)^{2k}\cdot h(z)$ with a holomorphic $h$ that does not vanish at $z_0$, so $\sqrt{h(z)}$ is holomorphic in a neighbourhood of $z_0$ and $\sqrt{f(z)} = (z-z_0)^k\sqrt{h(z)}$ is a well-defined meromorphic function in a neighbourhood of $z_0$), so only the zeros and poles of odd multiplicity are actual branch-points. [The situation is similar for $k$-th roots, zeros and poles whose multiplicity is a multiple of $k$ aren't branch-points, zeros and poles of other multiplicities are; essential singularities of $f$ may or may not be branch-points of $\sqrt[k]{f(z)}$, consider $\exp \frac{1}{z}$ and $z\exp \frac{1}{z}$.] Then you need to connect each branch-point to at least one other branch-point with a branch-cut to obtain a domain on which a holomorphic branch of $\sqrt[k]{f(z)}$ is defined.
In the situation here, $\infty$ is a pole of multiplicity $2$ of $z^2-1$, so we have only the two branch-points $1$ and $-1$ of $\sqrt{z^2-1}$. Any curve connecting the two branch-points gives a region on which a holomorphic branch of $\sqrt{z^2-1}$ is defined, whether that curve passes through $\infty$ or not. The most common choices of branch-cuts are the interval $[-1,1]$ and the other arc of the real axis (plus $\infty$) connecting the two points, $\bigl(\mathbb{R}\setminus (-1,1)\bigr) \cup \{\infty\}$ [the latter branch-cut makes $\sqrt{z^2-1}$ an even function on its domain].
The image shows the branch cuts in red for the first identity on the left and the second on the right.
I don't think splitting the radicand into factors is a good idea here. For $\sqrt{z+1}$ and $\sqrt{z-1}$ (resp. $\sqrt{1+z}$ and $\sqrt{1-z}$), you need a branch-cut connecting $\mp 1$ to $\infty$, but for $\sqrt{z^2-1}$, $\infty$ is in no way a special point. The branch-cuts of a product can [often] be chosen quite differently from the branch-cuts of the factors.
Why is there a branch cut extending the imaginary axis in the second case?
There shouldn't be. It's most likely an artifact of how the values are computed by the software.
How would I tailor the inbuilt python function to choose my own branch cuts?
Write a wrapper that multiplies the result of the inbuilt square root function with $-1$ for $z$ in the appropriate parts of the plane, according to your chosen branch-cuts.
Best Answer
In THIS ANSWER, I discussed the meaning of the identity
$$\log(z_1z_2)=\log(z_1)+\log(z_2) \tag 1$$
In that expression, the equality is interpreted as a set equality. This means for any value of $\log(z_1z_2)$ can be expressed as the sum of some value of $\log(z_1)$ and some value of $\log(z_2)$. In addition, the sum of any values of $\log(z_1)$ and $\log(z_2)$ can be expressed as some value of $\log(z_1z_2)$.
Now, suppose $f(z)=\sqrt{z^2-1}=\sqrt{(z-1)(z+1)}$. By definition, we have
$$\begin{align} f(z)&=e^{\frac12\log((z-1)(z+1))}\\\\ &=e^{\frac12\left(\log(z-1)+\log(z+1)\right)}\\\\ &=e^{\frac12\sqrt{z-1}}e^{\frac12\sqrt{z+1}}\\\\ \sqrt{z^2-1}&=\sqrt{z-1}\sqrt{z+1}\tag2 \end{align}$$
The equality in $(2)$ is a set equivalence analogous with $(1)$.
EXAMPLE:
For the example given in the OP, $z=-2$. We denote by $z_1$ and $z_2$, $z_1=z+1$ and $z_2=z-1$. Clearly, $z_1=-1$, $z_2=-3$, and $z_1z_2=3$.
The multi-valued term $\log(z_1z_2)$ is given by
$$\log(z_1z_2)=\log(3)=\log(|3|)+i2n\pi\tag 3$$
for any integer $n$. If we define $\log(z_1)=i\pi$ and $\log(z_2)=\log(|3|)+i\pi$, and if $n=1$ in $(3)$, then $\log(z+1)+\log(z-1)=\log(z^2-1)$. However, for any other $n$, the equality does not hold.
If we use $(3)$ to calculate $\sqrt{z^2-1}=\sqrt{z_1z_2}$, then we obtain
$$\sqrt{z_1z_2}=\sqrt{3}=\sqrt{|3|}e^{in\pi}\tag 4$$
For $n$ odd, the equality $\sqrt{z_1z_3}=-\sqrt{|3|}=\sqrt{z_1}\sqrt{z_2}$ holds, while for $n$ even, the equality does not hold.
Now, let's be more general. We cut the plane using the principal branch for $\sqrt{z-1}$ and observe that
Similarly, using the principal branch for $\sqrt{z+1}$, we see that
If we wish to use the principal branch for $\sqrt{z^2-1}$, then the equality $\sqrt{z^2-1}=\sqrt{z-1}\sqrt{z+1}$ does not hold in general. Rather, the equality holds as follows:
More simply expressed, the equality holds by using the expression $\arg(z^2-1)=\text{Arg}(z-1)+\text{Arg}(z+1)$.
Let's use the previous example in which $z=-2$. We find that $\text{Arg}(z-1)=\text{Arg}(z+1)=\pi$ and $\text{Arg}(z^2-1)=0$. Indeed, we have
$$\sqrt{z-1}\sqrt{z+1}=\sqrt{z^2-1}$$
when we take $\arg(z^2-1)=\text{Arg}(z^2-1)+2\pi=0+2\pi=2\pi$.
More simply, we have equality with $\arg(z^2-1)=\text{Arg}(z-1)+\text{Arg}(z+1)=2\pi$.