Or, "why is it possible to use trig sub to find the electric field along a uniformly charged rod?"
When I was taking integral calculus, the professor explained that one of the motivations for trigonometric substitution was to calculate the electric field at a point along the axis of a uniform line of charge. Such a calculation could be expressed by some variant of:
$$\int_{}^{}\frac{1}{\sqrt{x^{2} + a^{2}}}dx$$
which could be solved by trig substitution.
Why is this the case? All I know is that the classic trig sub integral comes from the equation of an ellipse:
$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$
But I don't know the link between the electric field along a line of charge and the ellipse equation, nor do I know why calculating the field results in a sum.
Okay, I can give you an answer, but I don't know that you'll like it. We'll start with Coulomb's law:
$$\vec E = \frac{q\hat r}{4\pi\epsilon_0r^2}$$
(This itself is a consequence of Gauss's law of electromagnetism, $\iint \vec E\cdot d\vec s = \iiint\frac{q}{\epsilon_0}dV$, considering a sphere with a point charge at the center and using the rotational symmetry to simplify.)
So in the case of the rod, you need to integrate,
$$\vec E = \int_{0}^{L}\frac{\lambda\hat r}{4\pi\epsilon_0r^2} dt = \int_{0}^{L}\frac{\lambda\vec r}{4\pi\epsilon_0r^3} dt$$
Note the parameters of integration; we're integrating now along the rod (which we're assuming, as the formula you gave seems to, goes along the x-axis with one end at the origin). So now sub in $\vec r = (a - t)\hat x + b\hat y$ to get
$$\vec E = \int_{0}^{L}\frac{\lambda((a-t)\hat x+b\hat y)}{4\pi\epsilon_0((a-t)^2 + b^2)^{\frac{3}{2}}}dt$$
Use substitution to get
$$\vec E = \int_{-a}^{L-a}\frac{\lambda(-t\hat x+b\hat y)}{4\pi\epsilon_0(t^2 + b^2)^{\frac{3}{2}}}dt$$
Separate the components and use another substitution to get
$$\hat y\int_{-a}^{L-a}\frac{\lambda b}{4\pi\epsilon_0(t^2 + b^2)^{\frac{3}{2}}}dt - \hat x\int_{-a}^{L-a}\frac{\lambda t}{4\pi\epsilon_0(t^2 + b^2)^{\frac{3}{2}}}dt = \hat y\int_{-a}^{L-a}\frac{\lambda b}{4\pi\epsilon_0(t^2 + b^2)^{\frac{3}{2}}}dt - \hat x\int_{a^2 + b^2}^{(L-a)^2 + b^2}\frac{\lambda}{8\pi\epsilon_0t^{\frac{3}{2}}}dt = \hat y\int_{-a}^{L-a}\frac{\lambda b}{4\pi\epsilon_0(t^2 + b^2)^{\frac{3}{2}}}dt + \hat x\frac{\lambda}{4\pi\epsilon_0}\left(\frac{1}{\sqrt{(L-a)^2 + b^2}} - \frac{1}{\sqrt{a^2 + b^2}}\right)$$
I don't think I'm wrong about the x-component, so I can only think that the formula you have there assumes that $a >> L$. But that's how you get there.
As for the resemblance to elliptic integrals, to the extent that question has an answer, it's because they both involve stretching Euclidean distances - the arc lengths of a circle in one case, the position of a point charge in the other.