Setup
Another possible approach is to consider the arc length parameterization of the two curves. This is a parameterization $\gamma : I \to \mathbb{R}^2$ (where $I$ is an interval) such that the distance a particle travels along the curve in $\mathbb{R}^2$ is equal to the amount of time it has traveled along the curve, i.e. the length of the curve from $\gamma(a)$ to $\gamma(b)$ is $b-a$.
An arc length parameterization of the circle, defined on the interval $[0,2\pi]$ is given by
$$ \gamma_C(t) = (\cos(t), \sin(t)).$$
Note that this parameterization can be extended by periodicity to the entire real line, but we don't really want to do that. For the square, an arc length parameterization on $[0,8]$ (where $8$ is the perimeter of the square) is given by
$$ \gamma_S(t) = \begin{cases}
(1,t) & \text{if $t \in [0,1]$,} \\
(2-t, 1) & \text{if $t \in (1,3]$,} \\
(-1,4-t) & \text{if $t \in (3,5]$,} \\
(t-6,-1) & \text{if $t \in (5,7]$, and} \\
(1,t-8) & \text{if $t \in (7,8]$.} \\
\end{cases}$$
Again, this can be extended by periodicity, but we don't need that here.
From the Square to the Circle
Observe that both of these functions are injective, which means that they are invertible. So a reasonable approach to this problem is to take a point on the square, map it back to the interval $[0,8]$, scale that interval down to $[0,2\pi]$, then map that to the circle. That looks like
$$ \gamma_C \circ \sigma_{\pi/4} \circ \gamma_S^{-1}, $$
where $\sigma_{a}$ denotes the scaling map
$$ \sigma_{a}(t) = at. $$
The inverse of $\gamma_S$ can be given explicitly by
$$\gamma_{S}^{-1}(x,y) = \begin{cases}
y & \text{if $x=1$ and $y\ge 0$,} \\
2-x & \text{if $x \in (-1,1)$ and $y=1$,} \\
4-y & \text{if $x=-1$ and $y\in [-1,1]$,} \\
6+x & \text{if $x \in (-1,1)$ and $y=-1$, and} \\
8+y & \text{if $x = 1$ and $y \in [-1,0)$.}
\end{cases}$$
As ugly as this is, this particular version of the mapping has the advantage of being relatively straightforward to implement in code.
From the Circle to the Square
The idea is similar: compute
$$ \gamma_S \circ \sigma_{4/\pi} \circ \gamma_C^{-1}. $$
The happy thing about the given arc length parameterization of the circle is that the arc length is equal to angle (in radians). So it is sufficient to determine the angle from the origin to a given point on the circle. Note that the slope of a line from the origin to a point $(x,y)$ is $\tan(\theta)$, where $\theta$ is the angle that line makes to the $x$-axis. This implies that
$$ \tan(\theta) = \frac{y}{x}. $$ It is tempting to then conclude that $$ \theta = \arctan\!\left( \frac{y}{x} \right), $$ but one needs to be a little careful about which quadrant things are living, and on which branch of the tangent function we want things to land. After a bit of work,
$$\gamma_{C}^{-1}(x,y) = \begin{cases}
\arctan(y/x) & \text{if $x \ge 0$ and $y \ge 0$, }\\
\arctan(y/x) + \pi/2 & \text{if $x < 0$, and} \\
\arctan(y/x) + 3\pi/4 & \text{if $x \ge 0$ and $y < 0$.}
\end{cases}$$
Alternatively, many programming languages also include a two-argument arctangent function (atan2), which returns an angle between $(-\pi,\pi]$. If one uses this function, $\gamma_C^{-1}$ could be reduced to
$$\gamma_C^{-1}(x,y) = \begin{cases}
\operatorname{atan2}(y,x) & \text{if $y \ge 0$, and} \\
\operatorname{atan2}(y,x) + 2\pi & \text{if $y < 0$.}
\end{cases}$$
