We leverage data-driven model discovery methods to determine the governing equations for the emergent behavior of heterogeneous networked dynamical systems. Specifically, we consider networks of coupled nonlinear oscillators whose collective behaviour approaches a limit cycle. Stable limit-cycles are of interest in many biological applications as they model self-sustained oscillations (e.g. heart beats, chemical oscillations, neurons firing, circadian rhythm). For systems that display relaxation oscillations, our method automatically detects boundary (time) layer structures in the dynamics, fitting inner and outer solutions and matching them in a data-driven manner. We demonstrate the method on well-studied systems: the Rayleigh Oscillator and the Van der Pol Oscillator. We then apply the mathematical framework to discovering lowdimensional dynamics in networks of semi-synchronized Kuramoto, Rayleigh, Rossler, and Fitzhugh-Nagumo oscillators, as well as heterogeneous combinations thereof. We also provide a numerical exploration of the dimension of collective network dynamics as a function of several network parameters, showing that the discovery of coarse-grained variables and dynamics can be accomplished with the proposed architecture.