The cortical surface atlases constructed from a large representative population of neonates are highly needed in the neonatal neuroimaging studies. However, existing neonatal cortical surface atlases are typically constructed from small datasets, e.g., tens of subjects, which are inherently biased and thus are not representative to the neonatal population. In this paper, we construct neonatal cortical surface atlases based on a large-scale dataset with 764 subjects. To better characterize the dynamic cortical development during the first postnatal weeks, instead of constructing just a single atlas, we construct a set of spatiotemporal atlases at each week from 39 to 44 gestational weeks. The central idea is that, for all cortical surfaces, we first group-wisely register them into the common space to ensure the unbiasedness. Then, rather than simply averaging over the co-registered cortical surfaces, which generally leads to over-smoothed cortical folding patterns, we adopt a spherical patch-based sparse representation using an augmented dictionary to overcome the noises and potential registration errors. Through the group-wise sparsity constraint, we obtain consistent geometric cortical folding attributes on the atlases. Our atlases preserve the sharp cortical folding patterns, thus leading to better registration accuracy when aligning new subjects onto the atlases.