I understand the need for phase unwrapping and the complexity behind it. But, I'm struggling to understand why the phase information gets wrapped between $-\pi$ to $\pi$ in the first place?
I'm not sure if I understand your problem, but I'll give it a try. If you have a complex frequency response
where $\phi(\omega)$ is the phase, then you probably know that for a given frequency $\omega_0$, the frequency response might as well be written as
where $k$ is some integer. The question is, when you compute the complex value $H(\omega_0)$, which of the infinitely many possible phases do you choose? The convention is to choose the principal value which is in the interval $(-\pi,\pi]$. The consequence is that whenever you compute the function $\phi(\omega)$ for continuous $\omega$, you'll get jumps as soon as the phase hits $-\pi$ or $\pi$, but no matter how you define the phase, you will not be able to avoid jumps. So if you want a continuous function you need to do phase unwrapping.
The math.h library atan2() function only has 2 scalar inputs, not the prior history required to know which way the input vector was rotating (assuming the input was a vector and a continuous function of something else, such as time). Without history information, a small turn in one direction looks exactly the same as a big turn in the other. So atan2() wraps. And computer programs use atan2() to calculate angles and phase from 2D vectors.