This paper presents a novel computational electromagnetics (CEM) technique, which hybridizes the periodic finite element method (FEM) with the method of moments (MoM), for efficient numerical modeling of electromagnetic scattering from metasurfaces consisting of truncated periodic or locally varying quasi-periodic array of structures. Based on the quasi-periodic nature of metasurfaces, the periodic FEM is employed to generate high-level macro basis functions (MBFs). Following that, a reduced MoM matrix is formed by using the MBFs with unknown coefficients belonging to each cell of the periodic array. The proposed hybrid method reduces the computational load significantly when compared with conventional numerical methods, especially for electrically large truncated metasurfaces involving arbitrarily inhomogeneous unit cells. Various numerical results are presented, and the accuracy and performance of the proposed method are tested against commercial codes, as well as an in-house FEM code.