add cubed-sphere

doc/background/spherical_samplings.rst

@@ -173,6 +173,19 @@ SHT.
 .. _geodesic grid: https://en.wikipedia.org/wiki/Geodesic_grid
 .. _icosahedral symmetry: https://en.wikipedia.org/wiki/Icosahedral_symmetry
 
+Subdivision of the cube (:class:`~pygsp.graphs.SphereCubed`)
+------------------------------------------------------------
+
+The sampling is made of a cube whose faces are subdivided into finer
+quadrilaterals to the desired resolution.
+The result is a convex polyhedron made of :math:`n=6⋅m^2` faces.
+The graph vertices represent the quadrilateral faces.
+
+The sampling is used in weather and climate modeling because edges and faces
+are of approximately equal length and area (better when
+``spacing='equiangular'``). It is hierarchical (for subdivisions arranged in
+powers of 2) but doesn't have a sampling theorem nor a fast SHT.
+
 Subdivision of the rhombic dodecahedron (:class:`~pygsp.graphs.SphereHealpix`)
 ------------------------------------------------------------------------------
 

pygsp/graphs/__init__.py

@@ -179,6 +179,7 @@
     SphereEquiangular
     SphereGaussLegendre
     SphereIcosahedral
+    SphereCubed
     SphereHealpix
     SphereRandom
 
@@ -219,6 +220,7 @@
     'SphereGaussLegendre',
     'SphereHealpix',
     'SphereIcosahedral',
+    'SphereCubed',
     'TwoMoons'
 ]
 

pygsp/graphs/nngraphs/spherecubed.py

@@ -0,0 +1,131 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+
+from pygsp.graphs import NNGraph  # prevent circular import in Python < 3.5
+from pygsp import utils
+
+
+class SphereCubed(NNGraph):
+    r"""Sphere sampled as a subdivided cube.
+
+    Background information is found at :doc:`/background/spherical_samplings`.
+
+    Parameters
+    ----------
+    subdivisions : int
+        Number of edges the cube's edges are divided into, resulting in
+        ``6*subdivisions**2`` vertices.
+    spacing : {'equiangular', 'equidistant'}
+        Whether the vertices (on a cube's face) are spaced by equal angles (as
+        in [4]_) or equal distances (as in [1]_).
+    kwargs : dict
+        Additional keyword parameters are passed to :class:`NNGraph`.
+
+    Attributes
+    ----------
+    signals : dict
+        Vertex position as latitude ``'lat'`` in [-π/2,π/2] and longitude
+        ``'lon'`` in [0,2π[.
+
+    See Also
+    --------
+    SphereEquiangular, SphereGaussLegendre : based on quadrature theorems
+    SphereHealpix, SphereIcosahedral : based on subdivided polyhedra
+    SphereRandom : random uniform sampling
+
+    Notes
+    -----
+    Edge weights are computed by :class:`NNGraph`. Gaussian kernel widths have
+    however not been optimized for convolutions on the resulting graph to be
+    maximally equivariant to rotation [7]_.
+
+    References
+    ----------
+    .. [1] R. Sadourny, Conservative finite-difference approximations of the
+       primitive equations on quasi-uniform spherical grids, 1972.
+    .. [2] EM O'Neill, RE Laubscher, Extended studies of a quadrilateralized
+       spherical cube Earth data base, 1976.
+    .. [3] R. A. White, S. W. Stemwedel, The quadrilateralized spherical cube
+       and quad-tree for all sky data, 1992.
+    .. [4] C. Ronchi, R. Iacono, P. Paolucci, The “Cubed-Sphere:” a new method
+       for the solution of partial differential equations in spherical
+       geometry, 1995.
+    .. [5] W. M. Putmana, S.-J. Lin, Finite-volume transport on various
+       cubed-sphere grids, 2007.
+    .. [6] https://en.wikipedia.org/wiki/Quadrilateralized_spherical_cube
+    .. [7] M. Defferrard et al., DeepSphere: a graph-based spherical CNN, 2019.
+
+    Examples
+    --------
+    >>> import matplotlib.pyplot as plt
+    >>> G = graphs.SphereCubed()
+    >>> fig = plt.figure()
+    >>> ax1 = fig.add_subplot(131)
+    >>> ax2 = fig.add_subplot(132, projection='3d')
+    >>> ax3 = fig.add_subplot(133)
+    >>> _ = ax1.spy(G.W, markersize=1.5)
+    >>> _ = G.plot(ax=ax2)
+    >>> G.set_coordinates('sphere', dim=2)
+    >>> _ = G.plot(ax=ax3, indices=True)
+
+    Equiangular and equidistant spacings:
+
+    >>> import matplotlib.pyplot as plt
+    >>> fig, ax = plt.subplots()
+    >>> graph = graphs.SphereCubed(4, 'equiangular')
+    >>> graph.set_coordinates('sphere', dim=2)
+    >>> _ = graph.plot('C0', edges=False, ax=ax)
+    >>> graph = graphs.SphereCubed(4, 'equidistant')
+    >>> graph.set_coordinates('sphere', dim=2)
+    >>> _ = graph.plot('C1', edges=False, ax=ax)
+    >>> _ = ax.set_title('Equiangular and equidistant spacings')
+
+    """
+
+    def __init__(self, subdivisions=3, spacing='equiangular', **kwargs):
+
+        self.subdivisions = subdivisions
+        self.spacing = spacing
+
+        def linspace(interval):
+            x = np.linspace(-interval, interval, subdivisions, endpoint=False)
+            x += interval / subdivisions
+            return x
+
+        a = np.sqrt(3) / 3
+        if spacing == 'equidistant':
+            x = linspace(a)
+            y = linspace(a)
+        elif spacing == 'equiangular':
+            x = a * np.tan(linspace(np.pi/4))
+            y = a * np.tan(linspace(np.pi/4))
+
+        x, y = np.meshgrid(x, y)
+        x, y = x.flatten(), y.flatten()
+        z = a * np.ones_like(x)
+
+        coords = np.concatenate([
+            [-y, x, z],   # North face.
+            [z, x, y],    # Equatorial face centered at longitude 0.
+            [-x, z, y],   # Equatorial face centered at longitude 90°.
+            [-z, -x, y],  # Equatorial face centered at longitude 180°.
+            [x, -z, y],   # Equatorial face centered at longitude 270°.
+            [y, x, -z],   # South face.
+        ], axis=1).T
+
+        coords /= np.linalg.norm(coords, axis=1)[:, np.newaxis]
+
+        super().__init__(coords, **kwargs)
+
+        lat, lon = utils.xyz2latlon(*coords.T)
+        self.signals['lat'] = lat
+        self.signals['lon'] = lon
+
+    def _get_extra_repr(self):
+        attrs = {
+            'subdivisions': self.subdivisions,
+            'spacing': self.spacing,
+        }
+        attrs.update(super(SphereCubed, self)._get_extra_repr())
+        return attrs

pygsp/graphs/nngraphs/spheregausslegendre.py

@@ -36,7 +36,8 @@ class SphereGaussLegendre(NNGraph):
     See Also
     --------
     SphereEquiangular : based on quadrature theorems
-    SphereIcosahedral, SphereHealpix : based on subdivided polyhedra
+    SphereIcosahedral, SphereCubed, SphereHealpix :
+        based on subdivided polyhedra
     SphereRandom : random uniform sampling
 
     Notes

pygsp/graphs/nngraphs/spherehealpix.py

@@ -46,7 +46,7 @@ class SphereHealpix(NNGraph):
     See Also
     --------
     SphereEquiangular, SphereGaussLegendre : based on quadrature theorems
-    SphereIcosahedral : based on subdivided polyhedra
+    SphereIcosahedral, SphereCubed : based on subdivided polyhedra
     SphereRandom : random uniform sampling
 
     Notes

pygsp/graphs/nngraphs/sphereicosahedral.py

@@ -42,7 +42,7 @@ class SphereIcosahedral(NNGraph):
     See Also
     --------
     SphereEquiangular, SphereGaussLegendre : based on quadrature theorems
-    SphereHealpix : based on subdivided polyhedra
+    SphereCubed, SphereHealpix : based on subdivided polyhedra
     SphereRandom : random uniform sampling
 
     Notes
@@ -121,6 +121,7 @@ def normalize(vertices):
             # Projecting pushes points away from the 12 base vertices, which
             # may make the point density more uniform.
             # See "A Comparison of Popular Point Configurations on S^2".
+            # As the equiangular vs equidistant spacing on the subdivided cube.
             normalize(mesh.vertices)
 
         if not dual:

pygsp/graphs/nngraphs/sphererandom.py

@@ -29,7 +29,8 @@ class SphereRandom(NNGraph):
     See Also
     --------
     SphereEquiangular, SphereGaussLegendre : based on quadrature theorems
-    SphereIcosahedral, SphereHealpix : based on subdivided polyhedra
+    SphereIcosahedral, SphereCubed, SphereHealpix :
+        based on subdivided polyhedra
     CubeRandom : randomly sampled cube
 
     References

pygsp/graphs/sphereequiangular.py

@@ -38,7 +38,8 @@ class SphereEquiangular(Graph):
     See Also
     --------
     SphereGaussLegendre : based on quadrature theorems
-    SphereIcosahedral, SphereHealpix : based on subdivided polyhedra
+    SphereIcosahedral, SphereCubed, SphereHealpix :
+        based on subdivided polyhedra
     SphereRandom : random uniform sampling
 
     Notes

pygsp/tests/test_graphs.py

@@ -572,6 +572,19 @@ def test_sphere_icosahedral(self, subdivisions=3):
         self._test_sphere(graph)
         self.assertEqual(np.sum(graph.signals['lat'] == 0), 4*2**subdivisions)
 
+    def test_sphere_cubed(self, subdivisions=9):
+        for spacing in ['equiangular', 'equidistant']:
+            graph = graphs.SphereCubed(subdivisions, spacing)
+            self.assertEqual(graph.n_vertices, 6*subdivisions**2)
+            self._test_sphere(graph)
+            n_pixels_at_equator = sum(abs(graph.signals['lat']) < 1e-10)
+            self.assertEqual(n_pixels_at_equator, 4*subdivisions)
+        # Ordering of the faces.
+        graph = graphs.SphereCubed(1, k=2)
+        coords = [[0, 0, 1], [1, 0, 0], [0, 1, 0],
+                  [-1, 0, 0], [0, -1, 0], [0, 0, -1]]
+        np.testing.assert_allclose(graph.coords, coords)
+
     def test_sphere_healpix(self, subdivisions=2):
         nside = 2**subdivisions
         graph = graphs.SphereHealpix(subdivisions)