1 files changed, 55 insertions, 0 deletions
diff --git a/tests/test_projections.py b/tests/test_projections.py
new file mode 100644
index 0000000..f1a5f90
--- /dev/null
+++ b/tests/test_projections.py
@@ -0,0 +1,55 @@
+import unittest
+
+import numpy as np
+
+from mapgen import projections as PJ
+
+
+class ProjectionTest(unittest.TestCase):
+    def test_equal_earth_roundtrip(self):
+        lat = np.array([0.0, 30.0, -60.0, 89.0, 10.0])
+        lon = np.array([0.0, 120.0, -170.0, 45.0, -179.0])
+        x, y = PJ.equal_earth_forward(lat, lon)
+        la, lo = PJ.equal_earth_inverse(x, y)
+        np.testing.assert_allclose(la, lat, atol=1e-6)
+        np.testing.assert_allclose(lo, lon, atol=1e-6)
+
+    def test_mollweide_roundtrip(self):
+        lat = np.array([0.0, 45.0, -80.0, 89.9])
+        lon = np.array([0.0, -100.0, 170.0, 20.0])
+        x, y = PJ.mollweide_forward(lat, lon)
+        la, lo = PJ.mollweide_inverse(x, y)
+        np.testing.assert_allclose(la, lat, atol=1e-6)
+        np.testing.assert_allclose(lo, lon, atol=1e-6)
+
+    def test_equal_area(self):
+        # equal lat/lon boxes at the equator and at 60° should differ in area by cos(lat) — as on the sphere
+        def box_area(f, lat0):
+            la = np.array([lat0, lat0, lat0 + 1, lat0 + 1]); lo = np.array([0.0, 1.0, 1.0, 0.0])
+            x, y = f(la, lo)
+            return 0.5 * abs(np.dot(x, np.roll(y, -1)) - np.dot(y, np.roll(x, -1)))
+        for f in (PJ.equal_earth_forward, PJ.mollweide_forward):
+            ratio = box_area(f, 60.0) / box_area(f, 0.0)
+            self.assertAlmostEqual(ratio, (np.sin(np.radians(61)) - np.sin(np.radians(60))) / np.sin(np.radians(1)),
+                                   delta=0.01)
+
+    def test_orthographic_centre_and_limb(self):
+        lat, lon, ok = PJ.orthographic_inverse(np.array([0.0, 0.99, 1.5]), np.array([0.0, 0.0, 0.0]), 40.0, 10.0)
+        self.assertAlmostEqual(lat[0], 40.0, places=9)
+        self.assertAlmostEqual(lon[0], 10.0, places=9)
+        self.assertTrue(ok[0] and ok[1] and not ok[2])
+
+    def test_reproject_image(self):
+        img = np.zeros((90, 180, 3), np.uint8)
+        img[:, :90] = (255, 0, 0)            # western hemisphere red
+        img[:, 90:] = (0, 0, 255)            # eastern blue
+        out = PJ.reproject(img, "equal_earth", 400)
+        h, w = out.shape[:2]
+        self.assertTrue(1.9 < w / h < 2.2)
+        self.assertEqual(tuple(out[h // 2, w // 4]), (255, 0, 0))
+        self.assertEqual(tuple(out[h // 2, 3 * w // 4]), (0, 0, 255))
+        self.assertEqual(tuple(out[2, 2]), PJ.BACKGROUND)      # outside the map outline
+        globe = PJ.globe(img, 0.0, -90.0, 200)
+        self.assertEqual(globe.shape[:2], (200, 200))
+        self.assertGreater(int(globe[100, 100, 0]), 150)        # centred on the red hemisphere
+        self.assertEqual(tuple(globe[2, 2]), PJ.BACKGROUND)