Multivariate random variables with scipy.stats rvs() function

The scipy.stats suite of statistical distributions (scipy.stats.norm, scipy.stats.uniform, scipy.stats.t etc) all produce univariate data series using their own .rvs() function, and only one has a multivariate rendition: multivariate_normal, which corresponds to numpy's numpy.random.randn((N,K)). In fact, virtually all of the statistical distributions found in numpy.random can produce multivariate data.

How can I extend the univariate distribution functions found in scipy.stats to multivariate number generation, given that it possesses some distributions not found in numpy.random like johnsonsu? Must I manually make a loop function myself that concatenates multiple univariates together? what should that look like