Abstract: Following an extensive investigation of various monolayer transition metal dichalcogenides (MX2), research interest has expanded to include multilayer systems. In bilayer MX2, the stacking order strongly impacts the local band structure as it dictates the local confinement and symmetry. Determination of stacking order in multilayer MX(2)domains usually relies on prior knowledge of in-plane orientations of constituent layers. This is only feasible in case of growth resulting in well-defined triangular domains and not useful in-case of closed layers with hexagonal or irregularly shaped islands. Stacking order can be discerned in the reciprocal space by measuring changes in diffraction peak intensities. Advances in detector technology allow fast acquisition of high-quality four-dimensional datasets which can later be processed to extract useful information such as thickness, orientation, twist and strain. Here, we use 4D scanning transmission electron microscopy combined with multislice diffraction simulations to unravel stacking order in epitaxially grown bilayer MoS2. Machine learning based data segmentation is employed to obtain useful statistics on grain orientation of monolayer and stacking in bilayer MoS2.