  We show that an optimized projection functions method can automatically
construct maximally localized Wannier functions even for bands with nontrivial
topology. We demonstrate this method on a tight-binding model of a
two-dimensional $\mathbb{Z}_2$ topological insulator, on a three-dimensional
strong $\mathbb{Z}_2$ topological insulator, as well as on first-principles
density functional theory calculated valence states of Bi$_2$Se$_3$. In all
cases, the resulting Wannier functions contain large imaginary components and
are more extended than those in the topologically trivial phase.
