worldgen

git clone https://git.godosa.eu/worldgen

master

history zip tar.gz

raw · 2638 bytes

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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)