In recent years twisted bi-layers of 2D materials became very popular in the field due to the possibility to totally change their electronic properties by simple rotation. At the same time, in the wide field of photonic crystals, this idea still remains almost untouched, and only some particular problems were considered. One of the reasons is the computational difficulty of the accurate consideration of Moiré superlattices that appear due to the superimposition of misaligned lattices. Indeed, the unit cell of the complex lattice is typically much larger than the original crystals and requires much more computational resources for the computations. Here, we propose a careful adaptation of the Fourier modal method in the form of the scattering matrices for the description of twisted 1D gratings' stacks. Our approach allows us to consider sublattices in close vicinity to each other and account for their interaction via the near-field. In the developed numerical scheme, we utilize the fact that each sublattice is only 1D-periodic and therefore simpler than the resulting 2D superlattice, as well as the fact that even a small gap between the lattices filters out high Fourier harmonics due to their evanescent origin. This accelerates the computations from 1 up to 3 and more orders of magnitude for typical structures depending on the number of harmonics. This paves the way for rigorous study of almost any photonic crystals of the proposed geometry and demonstration of specific Moiré-associated effects.