Generate Copula-correlated Samples With Specified Marginals In Python

I have N random variables (X1,...,XN) each of which is distributed over a specific marginal (normal, log-normal, Poisson...) and I want to generate a sample of p joint realizations of these variables Xi, given that the variables are correlated with a given Copula, using Python 3. I know that R is a better option but i want to do it in Python.

Following this method I managed to do so with a Gaussian Copula. Now I want to adapt the method to use a Archimedean Copula (Gumbel, Frank...) or a Student Copula. Ath the beginning of the Gaussian copula method, you draw a sample of p realizations from a multivariate normal distribution. To adapt this to another copula, for instance a bivariate Gumbel, my idea is to draw a sample from the joint distribution of a bivariate Gumbel, but I am not sure on how to implement this.

I have tried using several Python 3 packages : copulae, copula and copulas all provide the noption to fit a particular copula to a dataset but do not allow to draw a random sample from a given copula.

Can you provide some algorithmic insight on how to draw multivariate random samples from a given Copula with uniform marginals?

The following code implements the Clayton and AMH copulas. Page 4 of this paper shows how other kinds of copulas can be implemented.

import random
import math
import scipy.stats as st
def clayton(theta, n):
    v=random.gammavariate(1/theta,1)
    uf=[random.expovariate(1)/v for _ in range(n)]
    return [(k+1)**(-1.0/theta) for k in uf]

def amh(theta, n):
    # NOTE: Use SciPy RNG for convenience here
    v=st.geom(1-theta).rvs()
    uf=[random.expovariate(1)/v for _ in range(n)]
    return [(1-theta)/(math.exp(k)-theta) for k in uf]