What units does SphPolygon.area return?

[gerritholl]

The units returned by SphPolygon.area are not documented and not obvious.

Code Sample, a minimal, complete, and verifiable piece of code

In this MCVE, I construct a quadrangle near (0, 0) with an edge of 0.001 degrees. At such a small area, we can approximate the Earth as flat, so I can compare the polygon area calculated with pyresample against a multiplication of the lengths of the edges of the quadrangle. However, the results are very different.

from pyproj import Geod
from numpy import deg2rad, array
from pyresample.spherical import SphPolygon
from math import tau
x = 0.001
pol = SphPolygon(deg2rad(array([[0., 0.], [0., x], [x, x], [x, 0.]])))
g = Geod(ellps="sphere")
print(pol.area()*1.5e8/(2*tau))
print((g.inv(0, 0, x, 0)[2] * g.inv(0, 0, 0, x)[2]) / 1e6)

Problem description

I don't know the units of the result. The units are not documented at https://pyresample.readthedocs.io/en/latest/api/pyresample.html?highlight=SphPolygon#pyresample.spherical.SphPolygon.area. It's clearly not km². I thought it might be steradian, but multiplying by the surface of the Earth then dividing by 4π (or 2τ) doesn't yield a consistent result either.

Expected Output

I expect an output that I can interpret with documented units. The surface area of the indicated polygon should be in the order of 0.012 km² (Euclidian approximation), but I can't find an explicable conversion factor between the two numbers.

Actual Result, Traceback if applicable

0.003636094926640598
0.01236430006718607

Versions of Python, package at hand and relevant dependencies

Pyresample 1.22 on Python 3.9.7.

[djhoese]

Skimming the paper referenced in the docstring I found:

The user can choose the
units in which the polygon's area will be reported by specifying the radius of
the sphere in appropriate units. For example, if the radius is specified in kilo-
meters, the area will be returned in square kilometers.

[gerritholl]

I haven't specified any radius. Maybe pyresample uses some default?

[djhoese]

The default radius for SphPolygon's is 1. You can either specify a different spheroid radius to the _init__ method or based on the math at the end of the area method you could take the result and multiply by <your_radius> ** 2:

 return (sum(alpha) - (len(self.lon) - 2) * np.pi) * self.radius ** 2

[gerritholl]

I see, I can pass my own radius. Is that a supported feature? I tend to be hesitant to rely on undocumented API.

[gerritholl]

Passing radius works:

from pyproj import Geod
from numpy import deg2rad, array
from pyresample.spherical import SphPolygon
from math import tau
x = 0.001
pol = SphPolygon(deg2rad(array([[0., 0.], [0., x], [x, x], [x, 0.]])),
                 radius=6371)
g = Geod(ellps="sphere")
print(pol.area())
print((g.inv(0, 0, x, 0)[2] * g.inv(0, 0, 0, x)[2]) / 1e6)

results in

0.01236428559047198
0.01236430006718607

which is a match.

[gerritholl]

I'll make a small PR improving the documentation.

[gerritholl]

I also made a typo in my original MCVE, the surface area should have been 5.1e8 not 1.5e8... without that typo, I would have had a match in my initial post! That goes to show it's time for the weekend...