I know there are some research about explicit equations for affine models in $\mathbb{A}^2$ of many modular curves over $\mathbb{Q}$, for example of $X_i(N), X(N)$ (where $i = 0, 1, 2$) for small $N$.
These equations are very useful in order to study elliptic curves. (Indeed there are so many papers which use them in crucial point.)
But are there similar researches about explicit equations for the canonical morphism $X_1(N) \to X_0(N) : [E, P] \mapsto [E, \left< P \right>]$?
It is possible to find equations for $\pi : X_1(N) \to X_0(N)$ using work of Yifan Yang. In the article Defining equations of modular curves, he explains an algorithm to obtain equations for the modular curves $X_1(N)$ and $X_0(N)$. Let me explain how to use his results to get an equation for $\pi$.
Yang shows that there exist two modular units $F,G$ on $X_1(N)$ such that $F$ and $G$ are regular away from the cusp at infinity, and the orders of the poles of $F$ and $G$ at this cusp are relatively prime. Then $F$ and $G$ generate the function field of $X_1(N)$ and using the $q$-expansions of $F$ and $G$, it is easy linear algebra to compute an equation $P(F,G)=0$ of $X_1(N)$. If $F$ has a pole of order $m$ and $G$ has a pole of order $n$, then the equation has degree $n$ in $F$, and $m$ in $G$. The same can be done for $X_0(N)$, giving an equation $Q(X,Y)=0$.
Now the question is to express the modular units $X$ and $Y$ in terms of $F$ and $G$. Since Yang's construction uses explicit modular units, we know the $q$-expansions of $X,Y$ and $F,G$. Every rational function on $X_1(N)$ is of the form $\sum_{j=0}^{m-1} R_j(F) G^j$, where the $R_j$ are rational functions. In the case of $X$, if we know that the $R_j$ have degree $\leq d$, then we can simply use the $q$-expansions to find a linear relation of the form \begin{equation*} (*) \qquad \sum_{i=0}^d a_i F^i X = \sum_{i=0}^d \sum_{j=0}^{m-1} b_{i,j} F^i G^j. \end{equation*} with $a_i, b_{i,j} \in \mathbb{Q}$, and deduce $X$ in terms of $F$ and $G$ (and proceeding similarly for $Y$). As a variant, we may also search for a relation of the form $A(F,G)X=B(F,G)$.
Once the computer detects a relation $A(F,G)X=B(F,G)$ as above, then we should certify it. To this end, it suffices to show that the function $R=A(F,G)X-B(F,G)$ is regular at all cusps of $X_1(N)$ above $\infty \in X_0(N)$, and is zero at $\infty \in X_1(N)$. Indeed, the function $R$ is already regular away from these cusps. The cusps of $X_1(N)$ above $\infty \in X_0(N)$ are acted on simply transitively by the diamond automorphisms $\langle \gamma \rangle : X_1(N) \to X_1(N)$ with $\gamma \in \Gamma_0(N)/\pm \Gamma_1(N)$, the last group being isomorphic to $(\mathbb{Z}/N\mathbb{Z})^\times / \pm 1$. Checking regularity at these cusps can be done using the transformation formulas for modular units (Proposition 2 in Yang's article). An alternative way is to bound the orders of the poles of $X$ at these cusps. In fact these orders are equal since $X$ comes from $X_0(N)$, so we already know them. Then we just need to check that the order of vanishing of $R$ at $\infty \in X_1(N)$ is greater than the sum of the orders of the other possible poles. This can be done by computing the $q$-expansion to enough accuracy.
It would be nice to give an a priori bound on the degrees of the $R_j$ (that is, a bound on $d$). It is not difficult to find a polynomial $H$ in $F$ such that $HX$ is regular away from infinity. The idea is to cancel the poles of $X$ at the cusps $\langle \gamma \rangle \infty \in X_1(N) \backslash \{\infty\}$, by multiplying by suitable powers of $F-F(\langle \gamma \rangle \infty)$. Yet is not clear to me that the resulting function $HX$ is a polynomial in $F,G$, because the curve $P(F,G)=0$ may be singular, so that its ring of regular functions may be a strict subring in the normalisation $\mathcal{O}(X_1(N) \backslash \infty)$. However, searching for a relation as in $(*)$ should work in practice, by increasing the value of $d$ progressively if necessary.
What this answer doesn't address is whether the pairs $(F,G)$ and $(X,Y)$ can be chosen consistently, that is, in such a way that the resulting equation for $X_1(N) \to X_0(N)$ is as simple as possible. In fact, once $(F,G)$ has been found, there is a standard choice, namely $X=\sum_\gamma F | \gamma$ and $Y=\sum_\gamma G | \gamma$, where the average is over the diamond automorphisms. These are modular funtions on $X_0(N)$ and as Yang explains, the functions $X$ and $Y$ satisfy the assumptions of his theorem, so they generate the function field of $X_0(N)$ and we can compute the resulting model of $X_0(N)$. Then the method above gives an equation for the morphism $X_1(N) \to X_0(N)$. I don't know whether this choice of $(X,Y)$ produces simpler equations.
With the same technique, we can also express $F|\gamma$ and $G|\gamma$ in terms of $F$ and $G$ (recall that the $q$-expansions of $F|\gamma$ and $G|\gamma$ can be obtained using the transformation formulas for modular units). As a bonus, this gives equations for the diamond automorphisms $\langle \gamma \rangle : X_1(N) \to X_1(N)$. Note that if $(\mathbb{Z}/N\mathbb{Z})^\times / \pm 1$ is cyclic, generated by $\gamma_0$, then we only need to compute $F|\gamma_0$ and $G|\gamma_0$ in terms of $F$ and $G$.
I haven't done experiments with this method, but since Yang's functions and equations are pretty simple, we may hope the same for the equation of the morphism $X_1(N) \to X_0(N)$.
