The sigmoid range $(0,1)$ is technically open, because no input value maps to $0$ or $1$. You can get arbitrarily close to $0$ or $1$ but never equal them.
Taking an exponent followed by a log can cause overflow issues during computation, because exponents grow quickly. Sigmoid also gets close to zero quickly and this can underflow (potentially rounding to 0). However, the end result of the combined function is likely to be something within normal calculation range, because of the log (this is similar to - but even more extreme - multiplying and dividing by a very large number).
Although your quote suggests:
it is best to write the negative log-likelihood as a function of z
There is no analytical form where $E = -log(\sigma(z))$ can be re-written as a simple function of z. In theory, it can be calculated numerically by expanding terms carefully (beyond me, so won't show here), and some libraries may include this kind of expansion. In practice, to avoid numerical instability, many machine learning libraries will simply cap the value like this: $E = -log( max(\epsilon, \sigma(z)) )$ with $\epsilon$ a small number, maybe $10^{-15}$
A more common issue where the calculations are done accurately in neural networks is when considering the gradient of a loss function. If you use a sigmoid output layer alongside a binary cross-entropy cost function $E = -(y log(\sigma(z)) + (1-y)log(1-\sigma(z))) $, then some of the terms cancel out, and then the gradient contribution from training is trivially $\frac{\partial{E}}{\partial z} = y - \sigma(z)$ - there is no need in that case to calculate any log values despite them being in the loss function.
