An area element in spherical coordinates is given by
$$\mathrm d A=r^2\sin\phi \mathrm d\phi \mathrm d\theta,$$
following the convention shown here:
Place your sphere at the origin and say your projection is orthogonal to the $x,y$-plane, i.e. parallel to the $z$-axis. You now have a pixel. Let its bottom-left corner (as seen in the $x,y$-plane) have coordinate $(nd,md)$ with $d$ being the pixel width and $|n|,|m|<\frac{r}{d}$. You now project this pixel onto the sphere and want to know what the corresponding area on the sphere is. To do this, simply integrate the above mentioned area element over the limits in spherical coordinates corresponding to the change in $x$ and $y$, being $x_1=nd\rightarrow x_2=(n+1)d$ and $y_1=md\rightarrow y_2=(m+1)d$ respectively. However, we want to find these limits in terms of $\theta$ and $\phi$, i.e. $\theta_1,\theta_2$ and $\phi_1,\phi_2$, so that we can calculate the projected area of pixel $n,m$ by
$$A_{nm}=\int_{\phi_1}^{\phi_2}\int_{\theta_1}^{\theta_2}\mathrm d A.$$
We have that
\begin{align}
\theta&=\arctan\left(\frac{y}{x}\right) \\
\phi&=\arccos\left(\frac{z}{r}\right).
\end{align}
Now, the $\theta$-part is easy enough, as
\begin{align}
\theta_1&= \arctan\left(\frac{m}{n}\right)\\
\theta_2&= \arctan\left(\frac{m+1}{n+1}\right).
\end{align}
However, the $\phi$-part requires knowledge of $z$, but since we are always only projecting onto the surface of the sphere, we know that
$$x^2+y^2+z^2=r^2 \implies z=\sqrt{r^2-x^2-y^2},$$
(note that I've chosen the positive solution... If I were you, I would rotate my coordinate-system each time you're working with pixels lying outside of the first quadrant)
which gives
\begin{align}
z_1&=r\sqrt{1-\frac{d^2}{r^2}(n^2+m^2)} \\
z_2&=r\sqrt{1-\frac{d^2}{r^2}((n+1)^2+(m+1)^2)}.
\end{align}
But with this, we can calculate the remaining limits:
\begin{align}
\phi_1&=\arccos\left(\sqrt{1-\frac{d^2}{r^2}(n^2+m^2)}\right)\\
\phi_2&=\arccos\left(\sqrt{1-\frac{d^2}{r^2}\left((n+1)^2+(m+1)^2\right)}\right),
\end{align}
such that the integral evaluates to
\begin{align}
A_{nm}&=r^2(\theta_2-\theta_1)(\cos\phi_1-\cos\phi_2)\\
&= r^2\left(\arctan\left(\frac{m+1}{n+1}\right)-\arctan\left(\frac{m}{n}\right)\right)\left(\sqrt{1-\frac{d^2}{r^2}(n^2+m^2)}-\sqrt{1-\frac{d^2}{r^2}\left((n+1)^2+(m+1)^2\right)}\right)
\end{align}
In conclusion, this is a formula for the area of the pixel with $x,y$-coordinates $(nd,md)$ projected onto a sphere that has center at the origin and radius $r.$ Do this for every pixel whose projection is completely on the sphere (the formula does not work otherwise) and you can find your desired ratio.
Doing concrete calculations with this requires care w.r.t. principal values and such, but as I mentioned, I'd strongly recommend rotating your coordinate system when working with the pixels lying in other quadrants of the $x,y$-plane than the first.
Do leave a comment if you have any questions.
I will stick to your notations $\phi$ for longitude, $\theta$ for latitude, although one finds rather often the inverse convention.
Given: a center $M_0$ (defined by spherical coordinates $(\phi_0,\theta_0)$) and a radius $R$ ($0<R<\pi/2$, measured as an arc on the unit sphere).
The set of points $M$ of the sphere which are interior to the spherical disk with center $M_0$ and radius $R$ is given by the following dot product constraint:
$$\vec{OM}.\vec{OM_0}>\cos(R)$$
which is equivalent, using classical spherical coordinates, to:
$$\begin{pmatrix}\cos(\phi)\cos(\theta)\\ \sin(\phi)\cos(\theta)\\ \sin(\theta)\end{pmatrix} .\begin{pmatrix}\cos(\phi_0)\cos(\theta_0)\\ \sin(\phi_0)\cos(\theta_0)\\ \sin(\theta_0)\end{pmatrix}>\cos(R)\tag{1}$$
This constraint can be written under the form:
$$\cos(\theta_0)\cos(\theta)\cos(\phi-\phi_0)+\sin(\theta_0)\sin(\theta)>\cos(R)\tag{2}$$
which is an implicit equation in $(\phi,\theta)$ depending upon three parameters $(\phi_0,\theta_0,R)$ that can be visualized with this Geogebra animation (play with the sliders !):
https://www.geogebra.org/calculator/eanc7njp
The top sliders $f$ and $t$ refer to the coordinates $\phi_0$ and $\theta_0$ resp. of center $M_0$.
Here are two examples:
Fig. 1: An "ordinary" circle centered in (0,0) with radius $\pi/4$ rendered as a kind of ellipse.
Fig. 2: A "limit case" image of a circle belonging to northern hemisphere, tangent to the equator, (almost) passing throughout North Pole, explaining the almost linear segment ranging approximately from $-\pi/2$ to $\pi/2$.
Remark 1: constraint (1) has been given in the same form by @blamocur.
Remark 2: formula (2) could have been obtained directly by using the spherical law of cosines..
Remark 3: Explicit equations of the circle can be found here.
Best Answer
Well, I guess I'll take a stab at this. This is definitely a calculus problem. To take the area between two curves, you want to take the integral of the greater function minus the lesser function. For the yellow area, the greater function is $y=1$. The lesser function will take some manipulation. The formula for the circle is:
$(x-4)^2+(y-4)^2=16$
$(y-4)^2=16-(x-4)^2$
$y-4=-\sqrt{16-(x-4)^2}$
$y=4-\sqrt{16-(x-4)^2}$
So our integral is $\int^3_2[1-(4-\sqrt{16-(x-4)^2}]dx$=$\int^3_2-3dx+\int^3_2\sqrt{16-(x-4)^2}dx$. The first integral is $-3x$ evaluated from 2 to 3, or in other words, $-3(3)-[-3(2)]=-9+6=-3$.
For the second half of that integral, we'll use the info from Andreas's comment. We'll perform a change of variable
$u=x-4,du=dx$
$\int^3_2\sqrt{16-(x-4)^2}dx=\int^{-1}_{-2}\sqrt{16-u^2}du=\frac12[u\sqrt{16-u^2}+16sin^{-1}\frac u4]^{-1}_{-2}=\frac12[(x-4)\sqrt{16-(x-4)^2}+16sin^{-1}\frac{x-4}4]^3_2$
That solves the yellow area. For the other 2, you'll want to know where the 2 functions cross.
$(x-4)^2+(1-4)^2=16$
$(x-4)^2=7$
$x=4-\sqrt7$
For the green area, the 2 functions are the same, but it's evaluated from $4-\sqrt7$ to 2. For the orange area, the greater function is $y=2$, but the lesser function changes. It should be easy to see, though, that the right half is a rectangle. The left half is integrated from 1 to $4-\sqrt7$. Also, the first half of the integral has changed from $\int-3dx$ to $\int-2dx=-2x$. So the total orange area is $2x-\frac12[(x-4)\sqrt{16-(x-4)^2}+16sin^{-1}\frac{x-4}4]$ evaluated from 1 to $4-\sqrt7$ plus $1[2-(4-\sqrt7)]$, or $\sqrt7-2$.