Let, $\displaystyle f(z) = \sum\limits_{k=0}^{\infty} a_kz^k$ for $z \in \mathbb{D}$. Then defining $d_r := \operatorname{diam}f(r\mathbb{D})$ we see that $d_r/r$ is a non decreasing function in $r$ (this follows by applying the max-modulus principle to $\displaystyle \frac{f(z) - f(wz)}{z}$ on the disk $z \in r\mathbb{D}$ where, $|w| = 1$). Therefore, letting $r \to 0^+$ we see that $$2|a_1| = 2|f'(0)| = \limsup\limits_{r \to 0^+} \frac{d_r}{r} \le d_1. \tag{1}$$ Then the equality $2|a_1| = d_1$ implies $d_r/r = 2|a_1|$ for all $r \in [0,1)$.
Also, from Schwarz lemma we must have $$|f(z) - f(-z)| \le \frac{d_r}{r}|z|, \, \text{ for } z \in r\mathbb{D} \tag{2}$$ and in particular $2|f'(0)| = 2|a_1| \le d_r/r$ (letting, $|z| \to 0^+$). I.e., the equality $2|a_1| = d_r/r$ corresponds to equality in Schwarz lemma. Hence, we have $$f(z) - f(-z) = \frac{d_r}{r}z = 2a_1z \tag{3}$$ (wlog, we can assume that the unimodular constant is $1$ and in particular $a_1$ is a real number).
Now, consider the function $g(\theta) := |f(e^{i\theta}z) - f(-z)|^2$ where, we have fixed a $z \in \partial r\mathbb{D}$ (i.e., $|z| = r$). Then from $(2)$ we know that $g(\theta)$ is maximized when $\theta = 0$. In particular we must have $g'(0) = 0$.
Now, substituting from relation $(3)$ we note that \begin{align*}g'(\theta) &= \frac{d}{d\theta}\left|f(e^{i\theta}z) - f(z) + 2a_1z\right|^2 \\&= \frac{d}{d\theta} \left[\left(f(e^{i\theta}z) - f(z) + 2a_1z\right) \left(\overline{f(e^{i\theta}z) - f(z) + 2a_1z}\right) \right] \\&= 2 \Re \left[ ie^{i\theta}zf'(e^{i\theta}z)\left(\overline{f(e^{i\theta}z) - f(z) + 2a_1z}\right) \right]. \tag{4}\end{align*}
That is $\displaystyle g'(0) = -2|z|^2a_1\Im \left[f'(z)\right] = 0$ for all $|z| = r$, and hence $\Im \left[f'(z)\right] = 0$ for $|z| = r$ implies $f'(z) \equiv a_1$ in $\mathbb{D}$. That is $f(z) = a_0 + a_1z $ is a linear function. $\square$
An alternative approach: (inspired by Conrad's solution)
Let us denote $\displaystyle N(r) := \frac{1}{\pi r^2} \int_{r \mathbb{D}} |f'(z)|^2\,dx\,dy$ for $r \in [0,1]$.
Now, note that $\lim\limits_{r \to 0^+} N(r) = |f'(0)|^2 > 0$ (since, $f'(0) \neq 0$ otherwise it is trivial), i.e., since $f$ is locally injective near origin by the Area formula we have $$\frac{\text{Area}(f(\mathbb{rD}))}{\pi r^2} = N(r) = \frac{1}{\pi r^2}\int_{\mathbb D} |f'(z)|^2 \,dx\,dy = \sum_{k=1}^{\infty} k|a_k|^2r^{2k-2}$$ for all $r$ small enough. Therefore, $N(r)$ is strictly increasing for small $r > 0$ unless $a_k = 0$ for all $k \ge 2$, i.e., $f$ is linear.
Coming back to the problem if we assume $f$ is not linear then for small $r > 0$, $$|f'(0)|^2 = N(0) < N(r) = \frac{\text{Area}(f(\mathbb{rD}))}{\pi r^2} \le \frac{\pi d_r^2}{4\pi r^2} = |f'(0)|^2 \tag{5}$$ where, the second inequality in $(5)$ is due to the isodiametric inequality followed by the equality established in eqn $(1)$. Contradiction!
Hence, $f$ must be linear. $\square$