Mathematically, the probability density function ($\operatorname{PDF}$) for $Z$ is given by the integral:
$$\operatorname{PDF}(Z) = \int \delta\left(Z - f_a(X,Y)\right) \operatorname{PDF}(X,Y)\operatorname{d}X \operatorname{d}Y.$$ If the transformation $$\left[\begin{array}{c} X \\ Y \end{array}\right] \rightarrow \left[\begin{array}{c} f_a(X,Y) \\ g_a(X,Y) \end{array}\right] $$ is one to one (if not, break it down into regions where it is, and sum) then:
$$\operatorname{PDF}(Z) = \int \frac{\operatorname{PDF}(X(f_a, g_a),Y(f_a, g_a))}{\left|\frac{\partial f_a}{\partial X} \frac{\partial g_a}{\partial Y} - \frac{\partial f_a}{\partial Y} \frac{\partial g_a}{\partial X} \right|}\operatorname{d}g_a. $$ Some of the simpler choices for $g_a$ include $g_a= X$ or $g_a = Y$.
As for how I would do this, given the vague setup, I would set up a Metropolis-Hastings algorithm to produce a chain of $(Z, g_a)$ pairs from the bivariate density $$\frac{\operatorname{PDF}(X(f_a, g_a),Y(f_a, g_a))}{\left|\frac{\partial f_a}{\partial X} \frac{\partial g_a}{\partial Y} - \frac{\partial f_a}{\partial Y} \frac{\partial g_a}{\partial X} \right|},$$ and just drop the auxiliary variable, $g_a$, when returning results. You may need to implement a "sample every $N^{\mathrm{th}}$ step" limitation to reduce sample correlations. If the numbers don't have to be generated on command, but can be generated as a batch, then you can use the emcee MCMC sampling package.
scipy
(or anywhere else), though; you'll need to implement your own. $\endgroup$ – Christian Clason Sep 14 '16 at 7:41