I'll start with your second question first. Often, it is simpler to use $i/(n+1)$ (or $(i-1/2)/n$) instead of $i/n$ to avoid boundary issues. Let's say you've sorted your data so $Z_1,\ldots, Z_n$ refers to the smallest up to the biggest values in the data. Then, $Z_i$ is a pretty good estimate of the $(i/n)$th quantile of the distribution. But $i/n$ and $i/(n+1)$ are really close values and typically the true quantiles of $i/n$ and $i/(n+1)$ are also really close. Thus, we can argue that $Z_i$ is a pretty good estimate of the $(i/(n+1))$th quantile of the distribution. So, why does all this matter? Consider the biggest value in the data set, $Z_n$. Would this ever be a good estimate of quantile $n/n=1$? Of course not. For a Weibull distribution, quantile $1$ is always equal to $\infty$ yet $Z_n$ is always finite. That said, $Z_n$ is still a good approximation of quantile $n/(n+1)$.
As for your primary question, "why estimate by least squares?", you simply need to understand a simple fact of Weibull distributions. The Weibull CDF is given by
$$F(x) = 1-\exp(-(x/\delta)^c)$$
If we rearrange this expression, we get
$$\ln(x) = \dfrac 1 c \ln(-\ln(1-F(x))) + \ln(\delta)$$
Now, if we define $Y_i = \ln(Z_i)$ and $X_i=\ln(-\ln(1-[i/(n+1)]))$, then we would expect to get the relationship, $Y_i = a X_i + b$, where $a=1/c$ and $b=\ln(\delta)$. All you have to do then is run a least squares regression of $X$ vs $Y$ and convert the parameter estimates into $c$ and $\delta$.