The design complexity of photonic crystals and periodic gratings has been continuously increasing, driven by exploration of their unique physical phenomena and widespread applications. However, existing approaches for scattering modeling of periodic structures, typically relying on commercial software, often encounter challenges when adapting to complex configurations, especially in the context of near-field analysis and frequency responses near resonance. Therefore, the development of a more efficient and general scattering modeling approach to overcome these limitations has emerged as a fundamentally crucial and intricate task. In this paper, an efficient surface integration equation (SIE)-based method is developed to model the scattering properties of arbitrary-shaped 2D gratings with 1D periodicity. The SIE is solved with a Nystrom approach, which incorporates a local correction scheme and a Gaussian-Legendre quadrature rule to deal with singularity points with improved accuracy. The efficient evaluation of periodic Green's functions is achieved by combining an imaginary wavenumber extraction technique with an integral transformation approach. Furthermore, an over-determined matrix equation is constructed by testing the SIE with redundant observation points to mitigate potential internal resonance phenomena. The proposed approach is assessed through various numerical examples involving scatterers of different shapes and arrangements to demonstrate its accuracy and efficiency. The transmissivity spectra and surface field results, considering both normal and grazing incidence, are computed and compared against Comsol simulations and an independent wave-expansion-based analytical method. The method proposed is found to be superior in accuracy and efficiency to Comsol, a commercial software application, especially when complicated evanescent modes are excited. The proposed method serves as an efficient tool for modeling scattering from periodic gratings with arbitrary shapes.