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.