Calculating the solid angle of a 'polygon' in spherical coordinates

The magnitude of the solid angle in steradians is:

$$ \Omega=\frac{A}{r^2} $$

where $A$ is the area of ​​the spherical region and $r$ is the radius of the sphere.

In the case of an spherical polygon:

$$ A=\left[\left(\begin{aligned}\sum_{k=1}^n\alpha_k\end{aligned}\right)-(n-2)\pi\right]r^2 $$

where $\alpha_k$ is the k-th radian angle.

Therefore, knowing the coordinates of the vertices in a counterclockwise direction:

$$ V_k=(r\cos\phi_k\cos\theta_k,\,r\cos\phi_k\sin\theta_k,\,r\sin\phi_k), \quad\text{with}\;(\phi_k,\theta_k)\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right]\times[0,2\pi) $$

the angle between the normals to the planes on which two consecutive sides lie is:

$$ n_1\equiv\frac{V_k\times V_{k+1}}{||V_k\times V_{k+1}||},\quad n_2\equiv\frac{V_k\times V_{k-1}}{||V_k\times V_{k-1}||} \quad\Rightarrow\quad\beta_k=\arccos(n_1\cdot n_2) $$

which if the inequality is verified:

$$ (n_1\times n_2)\cdot V_k\ge 0 $$

we have $\alpha_k=\beta_k$, otherwise $\alpha_k=2\pi-\beta_k$.

The best option is to divide the region in triangles.

First of all we need to assume that the circular segment that join two points is part of a great circle.

That's the easiest assumption.

If the region is not convex, we can always add extra areas to be regarded as negative areas.

The area of a triangle on a unit sphere is simply (Girard's theorem): $$Area = A+B+C-\pi,$$ where $A,B,C$ are the angles between 2 points.

The angle between 2 points $p$ and $p'$ is given by $${\displaystyle {\begin{aligned}{A(p, p') }&={2 \arcsin \frac{1}{\sqrt{2}} \sqrt {1-\sin {\theta }\sin {\theta '}\cos {(\varphi -\varphi ')}-\cos {\theta }\cos {\theta '}}}\end{aligned}}}$$

For this last formula, check https://en.wikipedia.org/wiki/Spherical_coordinate_system#Distance_in_spherical_coordinates.

Once we have the distance between 2 points on a unit sphere, this is the chord on a great circle, a unit circle, and taking the angle is quite trivial: $2 \arcsin(D/2)$

I ended up going with this approach from the Wikipedia article as it is very simple to implement.

It is based on the area of the canonical spherical trapezoid you would imagine if you were applying trapezoidal integration on a sphere.

$\tan (0.5\,\Omega_i) = \tan ( 0.5\, (\theta_{i+1} - \theta_i))\, \sin( 0.5 (\phi_{i+1} + \phi_i))\, /\, \cos( 0.5\, (\phi_{i+1} - \phi_i))$

This is a signed area, so it will handle non-convex polygons automatically. It will also handle overlapping polygons with holes in them if the orientation of the polygons is the opposite for holes.