My attempt.
$f$ is an eigenform with integer coefficients $\in S_2(\Gamma_0(N))$, otherwise it is usually not true that its periods form a lattice in $\Bbb{C}$. As $f=\sum_{d|N} c_d \tilde{f}(dz)$ for some newform of lower level, we can restrict to the case that $f$ is a newform.
Following Cremona, taking finitely many $\alpha\in \Bbb{Q}$ representing the finitely many cusps, for some prime $p\equiv 1\bmod N$ not dividing the numerator and denominator of any $\alpha$, $$(1+p-a_p(f))\int_\alpha^{i\infty} f(z)dz=\int_\alpha^{i\infty} (1+p-T_p) f(z)dz$$ $$=\int_\alpha^{p\alpha} f(z)dz+\sum_{k=0}^{p-1} \int_\alpha^{(\alpha+k)/p} f(z)dz\tag{1}$$ $p\alpha$ and $(\alpha+k)/p$ are $\Gamma_0(N)$ equivalent to $\alpha$, whence the RHS is equal to a sum of integrals over closed-loops in $X_0(N)$.
From the convergence of the Rankin Selberg integral $\langle f E_2,f\rangle$ you'll get that $p+1-a_p(f)\ne 0$ for some $p\equiv 1\bmod N$ not dividing the numerator and denominator of the $\alpha$, obtaining that the subgroup of $\Bbb{C}$ generated by the $\int_\alpha^\beta f(z),\alpha,\beta\in \Bbb{Q}\cup i\infty$ is a lattice iff $\{\int_\gamma f(z),\gamma$ closed-loop in $X_0(N)\}$ is a lattice.
Then the theorem we need is that with $g=\dim_\Bbb{C} S_2(\Gamma_0(N))$ then $ \pi_1(X_0(N))^{ab}\cong \Bbb{Z}^{2g}$.
$S_2(\Gamma_0(N))$ has a $\Bbb{C}$-basis $f_1,\ldots,f_g$ of modular forms with integer coefficients. We choose it such that $f_1=f$. With $\gamma_1,\ldots,\gamma_{2g}$ a $\Bbb{Z}$-basis of $\pi_1(X_0(N))^{ab}$, and $\lambda_l=(\int_{\gamma_l} f_1(z)dz),\ldots,\int_{\gamma_l} f_g(z)dz))$, using that a harmonic function on $X_0(N)$ must be constant,you'll get that the columns of the $2g\times 2g$ matrix with rows $(\Re(\lambda_l),\Im(\lambda_l))$ are $\Bbb{R}$-linearly independent, whence so are its rows, so the $\lambda_l$ generate a lattice $\Lambda$ in $\Bbb{C}^g$.
Let $\Bbb{T}$ be the $\Bbb{Z}$-algebra generated by the Hecke operators. Using that $f_1$ is a newform we can take an element $P\in \Bbb{T}$ such that $Pf_j=0$ except for $Pf_1=cf_1\ne 0$.
Each $\int_\gamma Pf_j(z)dz$ is given by integrating $f_j$ on some curves from cusps to cusps, with the same $p$ as above and the second equality of $(1)$ you'll get that there is some $d$ such that for all closed-loop $\gamma$ in $X_0(N)$, $$(\int_\gamma (1+p-T_p)Pf_1(z)dz,\ldots,\int_\gamma (1+p-T_p)Pf_g(z)dz) \in \frac1d\Lambda$$
This proves that for any $\lambda\in \Lambda$, $((1+p-a_p(f))c\lambda_1,0\ldots,0) \in \frac1d\Lambda$, and hence $\{ \int_\gamma f_1(z)dz, \gamma \in \pi_1(X_0(N))\}$ is a lattice in $\Bbb{C}$, from which the $\int_\alpha^\beta f_1(z)dz,\alpha,\beta\in \Bbb{Q}\cup i\infty$ generate a lattice in $\Bbb{C}$.