First of all, the proper description of your map is as a holomorphic map ${\mathbb D}\to CP^1= {\mathbb C}\cup \{\infty\}$ (this is what it means to have poles) which is automorphic with respect to a certain group $\Gamma$ of linear fractional transformations of ${\mathbb D}$ preserving the order order 7 tiling:
$$
f(\gamma z)= f(z), \gamma\in \Gamma, z\in {\mathbb D}.
$$
Such $\Gamma$ is not unique and neither is $f$ (even after your normalization), unless you impose further restrictions. Firstly, you want to have $\Gamma$ which is the $(2,3,7)$ von Dyck group. To describe this group geometrically, pick one of the triangles $T$ in your tiling. Let $v$ be one of its vertices, $o$ the center of $T$ and $p$ the midpoint of an $e$ edge containing $v$. The group $\Gamma$ is discrete, generated by order 7 primitive elliptic rotation around $v$, order 2 elliptic rotation around $p$ and order $3$ rotation around $o$. (Rotations are understood here as linear fractional transformations conjugate to the euclidean rotations.) The quotient $X={\mathbb D}/\Gamma$ is topologically a 2-sphere with three distinguished points $[v], [o], [p]$, projections of $v, o, p$. The complex structure from ${\mathbb D}$ descends to a complex structure on $X$ making it a Riemann surface. As any compact Riemann surface of genus $0$, it is biholomorphic to the extended complex plane, $CP^1$. Such a biholomorphic map $h$ is not unique, but it becomes unique if you require
$$
h([o])=\infty, h([v])=0, h([p])=1.
$$
The composition of the (holomorphic) projection $\pi: {\mathbb D}\to X$ with $h$ is your automorphic function $f$. It will have a triple pole at $o$, order 7 zero at $v$ and also will send the edge $e$ to the unit segment $01$ in the real line. By automorphicity, the same will be true for all other vertices are centers of triangles.
This description is not explicit of course. However, the function $f$ (or its multivalued inverse) can be found as a solution of some ODE of the form $S(F)=\varphi$, where $\varphi$ is a certain meromorphic quadratic differential and $S$ is the Schwarzian derivative. Equivalently, $F$ is the ratio of two independent solutions of Fuchsian linear 2nd order ODE, of the type: $y'' + \varphi y=0$, where $y$ is a function of one complex variable.
This ODE can be made very explicit if you look at the target space (the extended complex plane); I can write it down if you are interested and when I have time. The differential $\varphi$ can be made very explicit on the target space but less so in the unit disk (one can find it as an infinite sum, a Poincare series). All of what I describe was known to Poincare and can be found in books on automorphic functions, like Ford's "Automorphic functions".
Edit. The inverse to the function $f$ is the ratio $F_1/F_2$, where $F_1(z), F_2(z)$ are two linearly independent solutions of the following holomorphic 2nd order ODE
$$
z(z-1) F'' + ((a+b+1)z-c)F' + abF=0
$$
on the complex plane, where
$$
|1-c|=1/7, |c-a-b|=1/2, |a-b|=1/3.
$$
One has to further fiddle with the choice of $F_1, F_2$ to make sure that $f^{-1}$ sends the upper half plane $H=\{z: Im(z)>0\}$ in the target complex plane to your favorite triangle $T$ (as above) in the unit disk.
See for instance Notes on differential equations and hypergeometric functions page 29. I do not have access to any math books right now and I got this reference just by googling "monodromy of hypergeometric functions".
Once you know $f: T\to H$, you get $f: {\mathbb D}\to CP^1$ by using the Schwarz reflection principle (to define $f$ on the image of $T$ under the hyperbolic reflection in the edge $e$) and then the automorphic equation $f(\gamma z)= f(z)$ on the rest of the unit disk.
