We leverage data-driven model discovery methods to determine governing equations for the emergent behavior of heterogeneous networked dynamical systems. Specifically, we consider networks of coupled nonlinear oscillators whose collective behavior approaches a limit cycle. Stable limit cycles are of interest in many biological applications, as they model self-sustained oscillations (e.g. heartbeats, 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 networks of semisynchronized Kuramoto, Rayleigh, Rossler, and FitzHugh-Nagumo oscillators, as well as heterogeneous combinations thereof, showing that the discovery of coarse-grained variables and dynamics can be accomplished with the proposed architecture.