This is probably very elementary, yet I am not sure how to approach this.
Suppose $f$ is a conformal map between Riemannian manifolds, which maps geodesics to geodesics. (i.e if $\alpha$ is a geodesic, then $f \circ \alpha$ is a geodesic).
Is it true that $f$ must be a "scaled isometry"? (Does there exist a constant $\lambda$ such that $df\in \lambda O(n)$ )
Note: We need to assume the domain is connected, of course.
It appears the answer is positive, i.e every conformal geodesics-preserving map is a scaled isometry. Here is a proof (based on Moishe Cohen's comment):
First, we prove this for the special case of the identity map $\operatorname{Id}:(M,g) \to (M,\tilde g)$. That is, suppose that every $g$-geodesic is also a $\tilde g$-geodesic, and that $\tilde g =e^{2f}g$ for some $f \in C^{\infty}(M)$.
Then,
$\tilde \nabla _X Y = \nabla _X Y + (X f )Y + (Y f )X - g(X,Y) \operatorname{grad}f \tag{1}$.
It's easy to see that the condition "$g$-geodesic $\Rightarrow$ $\tilde g$-geodesic" is equivalent to $\tilde \nabla _X X = \nabla _X X$ for every vector field $X$.
Combining this with $(1)$, we get $$2(X f )X - g(X,X) \operatorname{grad}f=0 \tag{2}$$ (for every $X$).
Now, let $p \in M$. We want to prove $\operatorname{grad}f(p) =0$.
Inserting $X=\operatorname{grad}f$ into $(2)$ we get
$$ \|\operatorname{grad}f\|^2 \operatorname{grad}f=0,$$
so $\operatorname{grad}f=0$. Since $M$ is connected $f$ must be constant.
Reduction of the general case to the case of the identity map:
Let $\phi:(M,g) \to (N,h)$ be a conformal geodesic map. First, suppose that $\phi$ is a diffeomorphism.
Since $\phi^{-1}:(N,h) \to (M,\phi^*h)$ is an isometry, it is confromal and geodesic. Hence, the composition
$$ \operatorname{Id}=\phi^{-1} \circ \phi:(M,g) \to (M,\phi^*h)$$ is confromal and geodesic, hence by the special case, it is a scaled isometry, i.e $\phi^*h=\lambda g$ for some $\lambda \in \mathbb{R}$. This shows $\phi$ is a scaled isometry.
If $\phi$ was not a diffeomorphism, we proceed as follows:
$\phi$ is a local diffeomorphism (inverse function theorem), hence it is a "local scaled isometry", in other words it is conformal with scaling factor which is a locally constant function (by the reasoning above for the diffeomorphism case). Since $M$ is connected, this scaling factor must be globally constant.