Add support for Galactic healpix files

README.md

@@ -76,7 +76,7 @@ data = fits.open('allsky.fits')[0].data
 def sampler(x, y):
     """
     x and y are arrays, giving the RA/Dec centers
-    for each pixel to extract
+    (in radians) for each pixel to extract
     """
     ... code to extract a tile from `data` here ...
 

toasty/tests/make_test_gal.py

@@ -0,0 +1,22 @@
+"""
+Resample test.hpx into Galactic coordinates, and create test_gal.hpx
+"""
+import numpy as np
+import healpy as hp
+from astropy.coordinates import Galactic, FK5
+import astropy.units as u
+
+nest = True
+map = hp.read_map('test.hpx', nest=nest)
+
+theta, phi = hp.pix2ang(64, np.arange(map.size), nest)
+l, b = phi, np.pi / 2 - theta
+
+g = Galactic(l, b, unit=(u.rad, u.rad))
+f = g.transform_to(FK5)
+ra, dec = f.ra.rad, f.dec.rad
+
+map = hp.get_interp_val(map, np.pi / 2 - dec, ra, nest)
+
+hp.write_map('test_gal.hpx', map, nest=nest, coord='G',
+             fits_IDL=False)

toasty/tests/test_toasty.py

@@ -29,7 +29,7 @@ def mock_sampler(x, y):
 def test_iter_tiles_path(depth):
     result = set(r[0] for r in iter_tiles(mock_sampler, depth))
     expected = set(['{n}/{y}/{y}_{x}.png'.format(n=n, x=x, y=y)
-                    for n in range(depth+1)
+                    for n in range(depth + 1)
                     for y in range(2 ** n)
                     for x in range(2 ** n)])
     assert result == expected
@@ -68,11 +68,11 @@ def test_reference_wtml():
     val = parseString(wtml)
 
     assert ref.getElementsByTagName('Folder')[0].getAttribute('Name') == \
-      val.getElementsByTagName('Folder')[0].getAttribute('Name')
+        val.getElementsByTagName('Folder')[0].getAttribute('Name')
 
     for n in ['Credits', 'CreditsUrl', 'ThumbnailUrl']:
         assert ref.getElementsByTagName(n)[0].childNodes[0].nodeValue == \
-          val.getElementsByTagName(n)[0].childNodes[0].nodeValue
+            val.getElementsByTagName(n)[0].childNodes[0].nodeValue
 
     ref = ref.getElementsByTagName('ImageSet')[0]
     val = val.getElementsByTagName('ImageSet')[0]
@@ -83,6 +83,7 @@ def test_reference_wtml():
 def cwd():
     return os.path.split(os.path.abspath(__file__))[0]
 
+
 def test_wwt_compare_sky():
     """Assert that the toast tiling looks similar to the WWT tiles"""
     direc = cwd()
@@ -113,11 +114,27 @@ def test_healpix_sampler():
         image_test(expected, result, "Failed for %s" % pth)
 
 
+@pytest.mark.skipif('not HAS_ASTRO')
+def test_healpix_sampler_galactic():
+
+    direc = cwd()
+
+    im = fits.open(os.path.join(direc, 'test_gal.hpx'))[1].data['I']
+
+    sampler = healpix_sampler(im, nest=True, coord='G')
+
+    for pth, result in iter_tiles(sampler, depth=1):
+        expected = read_png(os.path.join(direc, 'test_sky', pth))
+        expected = expected[:, :, 0]
+
+        image_test(expected, result, "Failed for %s" % pth)
+
+
 @pytest.mark.skipif('not HAS_ASTRO')
 def test_guess_healpix():
     pth = os.path.join(cwd(), 'test.hpx')
     d, nest, coord = tile._guess_healpix(pth)
-    assert nest == True
+    assert nest is True
     assert coord == 'C'
     np.testing.assert_array_equal(d, fits.open(pth)[1].data['I'])
 
@@ -138,8 +155,8 @@ def null_merge(mosaic):
             assert im.max() != 0
 
 
-
 class TestToaster(object):
+
     def setup_method(self, method):
         self.base = mkdtemp()
         self.cwd = cwd()
@@ -172,8 +189,6 @@ def test_no_merge(self):
         self.verify_toast()
 
 
-
-
 reference_wtml = """
 <Folder Name="ADS All Sky Survey">
 <ImageSet Generic="False" DataSetType="Sky" BandPass="Visible" Name="allSources_512" Url="allSources_512/{1}/{3}/{3}_{2}.png" BaseTileLevel="0" TileLevels="3" BaseDegreesPerTile="180" FileType=".png" BottomsUp="False" Projection="Toast" QuadTreeMap="" CenterX="0" CenterY="0" OffsetX="0" OffsetY="0" Rotation="0" Sparse="False" ElevationModel="False">

toasty/tile.py

@@ -13,15 +13,16 @@
 from collections import defaultdict, namedtuple
 
 level1 = [[np.radians(c) for c in row]
-          for row in  [[(0, -90), (90, 0), (0, 90), (180, 0)],
-                       [(90, 0), (0, -90), (0, 0), (0, 90)],
-                       [(0, 90), (0, 0), (0, -90), (270, 0)],
-                       [(180, 0), (0, 90), (270, 0), (0, -90)]]
-                       ]
+          for row in [[(0, -90), (90, 0), (0, 90), (180, 0)],
+                      [(90, 0), (0, -90), (0, 0), (0, 90)],
+                      [(0, 90), (0, 0), (0, -90), (270, 0)],
+                      [(180, 0), (0, 90), (270, 0), (0, -90)]]
+          ]
 
 Pos = namedtuple('Pos', 'n x y')
 Tile = namedtuple('Tile', 'pos increasing corners')
 
+
 def _postfix_corner(tile, depth, bottom_only):
     """
     Yield subtiles of a given tile, in postfix order
@@ -134,18 +135,17 @@ def iter_tiles(data_sampler, depth, merge=True):
     (pth, tile) : str, ndarray
       pth is the relative path where the tile image should be saved
     """
-    if merge == True:
+    if merge is True:
         merge = _default_merge
 
     parents = defaultdict(dict)
 
     for node, c, increasing in iter_corners(max(depth, 1),
-                                               bottom_only=merge):
+                                            bottom_only=merge):
 
         l, b = subsample(c[0], c[1], c[2], c[3], 256, increasing)
         img = data_sampler(l, b)
 
-
         for pth, img in _trickle_up(img, node, parents, merge, depth):
             yield pth, img
 
@@ -330,7 +330,7 @@ def _guess_healpix(pth, extension=None):
     return data, nest, coord
 
 
-def healpix_sampler(data, nest=False, coord = 'C', interpolation='nearest'):
+def healpix_sampler(data, nest=False, coord='C', interpolation='nearest'):
     """
     Build a sampler for Healpix images
 
@@ -354,6 +354,8 @@ def healpix_sampler(data, nest=False, coord = 'C', interpolation='nearest'):
     of (lon, lat)
     """
     from healpy import ang2pix, get_interp_val, npix2nside
+    from astropy.coordinates import Galactic, FK5
+    import astropy.units as u
 
     interp_opts = ['nearest', 'bilinear']
     if interpolation not in interp_opts:
@@ -362,14 +364,16 @@ def healpix_sampler(data, nest=False, coord = 'C', interpolation='nearest'):
     if coord.upper() not in 'CG':
         raise ValueError("Invalid coord %s. Must be 'C' or 'G'" % coord)
 
-    #XXX Add support for G
-    if coord.upper() == 'G':
-        raise NotImplementedError("coord='G' not yet supported")
-
+    galactic = coord.upper() == 'G'
     interp = interpolation == 'bilinear'
     nside = npix2nside(data.size)
 
     def vec2pix(l, b):
+        if galactic:
+            f = FK5(l, b, unit=(u.rad, u.rad))
+            g = f.transform_to(Galactic)
+            l, b = g.l.rad, g.b.rad
+
         theta = np.pi / 2 - b
         phi = l
 
@@ -403,7 +407,7 @@ def vec2pix(l, b):
         l[l < 0] += 2 * np.pi
         l = nx * (1 - l / (2 * np.pi))
         l = np.clip(l.astype(np.int), 0, nx - 1)
-        b = ny * (1 - (b + np.pi/2) / np.pi)
+        b = ny * (1 - (b + np.pi / 2) / np.pi)
         b = np.clip(b.astype(np.int), 0, ny - 1)
         return data[b, l]