add area() helper method

README.md

@@ -287,6 +287,7 @@ b'$EIGNQ,RMC*24\r\n'
  - `haversine` - finds great circle distance in km between two sets of (lat, lon) coordinates
  - `planar` - finds planar distance in m between two sets of (lat, lon) coordinates
  - `bearing` - finds bearing in degrees between two sets of (lat, lon) coordinates
+ - `area` - finds spherical area bounded by two sets of (lat, lon) coordinates
 
 See [Sphinx documentation](https://www.semuconsulting.com/pynmeagps/pynmeagps.html#module-pynmeagaps.nmeahelpers) for details.
 

RELEASE_NOTES.md

@@ -4,6 +4,7 @@
 
 CHANGES:
 
+1. Add area() helper method to calculate spherical area of bounding box.
 1. Sphinx documentation and docstrings enhanced to include global constants and decodes.
 1. `socket_stream.SocketStream` class renamed to `socket_wrapper.SocketWrapper` class for clarity.
 

examples/utilities.py

@@ -13,6 +13,7 @@
 from datums import DATUMS  # assumes this is in same folder
 
 from pynmeagps import (
+    area,
     bearing,
     ecef2llh,
     haversine,
@@ -44,21 +45,21 @@
     f" {LAT2}, {LON2} using default WGS84 datum...",
 )
 dist = haversine(LAT1, LON1, LAT2, LON2)
-print(f"Distance: {dist} km")
+print(f"Distance: {dist:.4f} km")
 
 print(
     f"\nFind planar distance between {LAT1}, {LON1} and",
     f" {LAT3}, {LON3} using default WGS84 datum...",
 )
 dist = planar(LAT1, LON1, LAT3, LON3)
-print(f"Distance: {dist} m")
+print(f"Distance: {dist:.4f} m")
 
 print(
     f"\nFind bearing between {LAT1}, {LON1} and",
     f" {LAT2}, {LON2} using default WGS84 datum...",
 )
 brng = bearing(LAT1, LON1, LAT2, LON2)
-print(f"Bearing: {brng} degrees")
+print(f"Bearing: {brng:.4f} degrees")
 
 X, Y, Z = 3822566.3113, -144427.5123, 5086857.1208
 print(f"\nConvert ECEF X: {X}, Y: {Y}, Z: {Z} to geodetic using default WGS84 datum...")
@@ -82,7 +83,14 @@
     f"{LAT2}, {LON2} using alternate {DATUM} ({ellipsoid}) datum...",
 )
 dist = haversine(LAT1, LON1, LAT2, LON2, a / 1000)
-print(f"Distance: {dist} km")
+print(f"Distance: {dist:.4f} km")
+
+print(
+    f"\nFind spherical area bound by {LAT1}, {LON1} and",
+    f"{LAT2}, {LON2}...",
+)
+bbarea = area(LAT1, LON1, LAT2, LON2)
+print(f"Area: {bbarea:.4f} km²")
 
 print(
     f"\nConvert ECEF X: {X}, Y: {Y}, Z: {Z} to",

src/pynmeagps/nmeahelpers.py

@@ -546,6 +546,29 @@ def bearing(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
     return brng
 
 
+def area(
+    lat1: float,
+    lon1: float,
+    lat2: float,
+    lon2: float,
+    radius: int = WGS84_SMAJ_AXIS / 1000,
+) -> float:
+    """
+    Calculate spherical area bounded by two coordinates.
+
+    :param float lat1: lat1
+    :param float lon1: lon1
+    :param float lat2: lat2
+    :param float lon2: lon2
+    :param float radius: radius in km (Earth = 6378.137 km)
+    :return: area in km²
+    :rtype: float
+    """
+
+    phi1, phi2 = [c * pi / 180 for c in (lat1, lat2)]
+    return pow(radius, 2) * pi * abs(sin(phi1) - sin(phi2)) * abs(lon1 - lon2) / 180
+
+
 def get_gpswnotow(dat: datetime) -> tuple:
     """
     Get GPS Week number (Wno) and Time of Week (Tow)

tests/test_static.py

@@ -14,6 +14,7 @@
 
 from pynmeagps import NMEAMessage, NMEAMessageError, NMEAReader, NMEATypeError
 from pynmeagps.nmeahelpers import (
+    area,
     bearing,
     calc_checksum,
     date2str,
@@ -488,6 +489,16 @@ def testbearing(self):
         res = bearing(-12.645, 34.867, -34.1745, 48.27846)
         self.assertAlmostEqual(res, 152.70835788275326, 4)
 
+    def testarea(self):
+        res = area(51.23, -2.41, 53.205, -2.34)
+        self.assertAlmostEqual(res, 1049.5657, 4)
+        res = area(53.48280729, -2.24225376, 53.46814647, -2.20192543)
+        self.assertAlmostEqual(res, 4.3606, 4)
+        res = area(53.69865772, -2.68269539, 53.22939103, -1.39218885)
+        self.assertAlmostEqual(res, 4467.6282, 4)
+        res = area(-12.645, 34.867, -34.1745, 48.27846)
+        self.assertAlmostEqual(res, 3264291.8230, 4)
+
     def testgpsweek(self):
         dats = [
             (2023, 1, 1),