Formula of this map (Escher's work)
Its formula is $$ f:D\to\mathbb{R}^3:(\phi,\theta)\mapsto(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta) $$ where $D=\{(\phi,\theta)\in[-\pi,\pi]\times[0,\pi]:\sin(4\phi+12\theta)>0\}$.
Here is a code in Mathematica
ParametricPlot3D[{Cos[\[Phi]] Sin[\[Theta]],
Sin[\[Theta]] Sin[\[Phi]],
Cos[\[Theta]]}, {\[Phi], -\[Pi], \[Pi]}, {\[Theta], 0, \[Pi]},
RegionFunction -> (Sin[4 #4 + 12 #5] > 0 &), Mesh -> None,
BoundaryStyle -> Black]
and the graph
It seems to me that the bands are supposed to twist indefinitely when they approach the poles (as observed by Erick Wong). That effect is achieved by using loxodromes as boundary lines. The loxodromes are roughly the spherical equivalents of the logarithmic spirals in that a loxodrome makes a constant angle with all the latitudes that it crosses - very much like the logarithmic spirals that have constant angles between their tangents and the radii.
The attaced image is generated by Mathematica. I used $m=1/2$ (see the link for the meaning of this parameter), so e.g. one of the bands was generated by
ParametricPlot3D[{Cos[t+u]/Cosh[t/2],Sin[t+u]/Cosh[t/2],Tanh[t/2]},{t,-8,8},{ u,0,Pi/4},PlotPoints->{81,5}]
and the other parts came out with the same formula, but the range of the variable $u$ shifted by an integer multiple of $\pi/2$. The range of $t$ is should be symmetric, but the endpoints (here $\pm8$) are, again, a largely arbitrary choice of mine. I did a linear change of parameters to that on the Wikipage in order to make sure that all the loxodromes cover the same interval of latitudes, and also to make one half of the boundaries of the tiles to have constant latitude.
Here's the image. To get a better match with Escher's painting it might be necessary to fine-tune the value of $m$ further.
To see the "endless twisting" here is a close-up of the polar region. The image is necessarily quite flat now. The formula is the same, but the range of the parameter $t$ is now $4\le t\le 16$.
Edit: Here's an animated version