The general solution of your equations in a simply connected domain on which $r_2\not=0$ and $r_1\not=\pm1$ is $$ \beta = \frac12 + \frac1{{(r_1}^2{-}1)}\, \left(\frac{\partial a}{\partial\theta_1}+b(\theta_1,r_2)\right) \quad\text{and}\quad \gamma= \frac12 + \frac1{{r_2}^2}\, \left(\frac{\partial a}{\partial\theta_2}+c(\theta_2,r_1)\right), $$ where $a = a(\theta_1,\theta_2)$ is a function of $\theta_1$ and $\theta_2$ only, $b$ is a function of $\theta_1$ and $r_2$ only, and $c$ is a function of $\theta_2$ and $r_1$ only.
To see this, set $\bar\beta = ({r_1}^2{-}1)\bigl(\beta-\tfrac12({r_1}^2{-}1)\bigr)$$\bar\beta = ({r_1}^2{-}1)\bigl(\beta-\tfrac12\bigr)$ and $\bar\gamma = {r_2}^2\bigl(\gamma-\tfrac12\bigr)$ and note that the given equations imply the constant coefficient linear equations $$ \frac{\partial\bar\beta}{\partial r_1} = \frac{\partial\bar\gamma}{\partial r_2} = \frac{\partial\bar\beta}{\partial\theta_2}-\frac{\partial\bar\gamma}{\partial\theta_1} = 0, $$ which are easily solved.
I do not see how you can choose $a$, $b$, and $c$ so that your 'boundary conditions' are satisfied. Because the general solution depends on three functions of two variables, one would expect to be able to specify (initial) conditions along a surface in the domain, but not along a hypersurface (which is what boundaries generally are).