  We present an all-electron, periodic {\GnWn} implementation within the
numerical atomic orbital (NAO) basis framework. A localized variant of the
resolution-of-the-identity (RI) approximation is employed to significantly
reduce the computational cost of evaluating and storing the two-electron
Coulomb repulsion integrals. We demonstrate that the error arising from
localized RI approximation can be reduced to an insignificant level by
enhancing the set of auxiliary basis functions, used to expand the products of
two single-particle NAOs. An efficient algorithm is introduced to deal with the
Coulomb singularity in the Brillouin zone sampling that is suitable for the NAO
framework. We perform systematic convergence tests and identify a set of
computational parameters, which can serve as the default choice for most
practical purposes. Benchmark calculations are carried out for a set of
prototypical semiconductors and insulators, and compared to independent
reference values obtained from an independent $G_0W_0$ implementation based on
linearized augmented plane waves (LAPW) plus high-energy localized orbitals
(HLOs) basis set, as well as experimental results. With a moderate (FHI-aims
\textit{tier} 2) NAO basis set, our $G_0W_0$ calculations produce band gaps
that typically lie in between the standard LAPW and the LAPW+HLO results.
Complementing \textit{tier} 2 with highly localized Slater-type orbitals
(STOs), we find that the obtained band gaps show an overall convergence towards
the LAPW+HLO results. The algorithms and techniques developed in this work pave
the way for efficient implementations of correlated methods within the NAO
framework.
