We must dispense the assumption that the square root of a positive integer is uniquely defined; indeed, $(-1)^2=1^2=1$, so that both $-1$ and 1 are valid roots of $1$.
Theoretically, this is an issue for every positive integer, but the reason that $\sqrt{1}=-1$ is important is that it offers the opportunity for $\frac{1+\sqrt{1+z^2}}{2}$ to be zero, at which point $\log \frac{1+\sqrt{1+z^2}}{2}$ will be undefined, and will become a branch point of $\log$. By defining $\sqrt{1+z^2}$ at $z=0$ to be $1$, we eliminate the possibility of there being a zero of $\frac{1+\sqrt{1+z^2}}{2}$, and thus $\infty,\pm i$ are the only remaining branch points.
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
With regard to 1) Calculus may shed a little light. If $f(z)=\sqrt{z}$, then $f'(z)=\dfrac{1}{2\sqrt{z}}$, and we see that $f'(z)$ is undefined at $z=0$ (where $f$ vanished). In terms of real variables the derivative doesn't exist at $0$, but in terms of complex variables $0$ is a non-isolated singularity, i.e. $f(z)$ is not analytic in a deleted neighborhood of $0$. The same holds true for $f(z)=\sqrt[n]{z}$.
It seems like 2) may indicate some misunderstanding related to an answer to 3). Geometrically a branch cut is a cut, a disconnection in the plane usually between pairs of branch points. It's as if a razor had been used. A richer geometric picture emerges from considering an alternative approach to branch cuts for handling branch points of multi-valued functions called Riemann Surfaces.
The encircling of branch points in the linked answer you referenced is a standard approach to determining branch cuts. The issue we must avoid is a discontinuity of the function after making a complete circuit around a point. We need to ensure that traversing a closed loop produces a change in the argument of $f(z)$ which is a multiple of $2\pi$. The simplest example is $f(z)=\sqrt{z}$. A closed loop around $z=0$ produces a change in the argument of $f(z)$ by $\pi$ for a change of $2\pi$ in the argument of $z$. That is what indicates that a branch cut from $0$ to $\infty$ is required. A closed loop around $0$ cannot be allowed, because it's impossible to avoid a discontinuity of $f(z)$. (Think about the loop $e^{it}$. At $t=0, z=1, \sqrt{1}=1\text{ and }\arg f(z) = 0$, but after a circuit around $z=0$ at $t=2\pi, z=1\text{ again, but, } \sqrt{1}=-1\text{ since}\arg f(z) = \pi$.
What the linked answer shows is that a loop around either $z=i$ or $-i$ alone produces a change of $\pi$ in $\arg \sqrt{z^2+1}$ which would be a problem, but a loop enclosing both of the points produces a change of $2 \pi$ in $\arg \sqrt{z^2+1}$. A branch cut joining $z=i$ and $z=-i$ ensures that any loop enclosing one of the points must enclose the other.