Added initial dask class to represent HiPS datasets

Author: astrofrog

reproject/hips/dask_array.py

@@ -0,0 +1,85 @@
+import os
+import urllib
+
+import numpy as np
+from astropy.io import fits
+from astropy.wcs import WCS
+from astropy_healpix import HEALPix, level_to_nside
+
+from .utils import is_url, load_properties, tile_filename
+
+
+class HiPSArray:
+
+    def __init__(self, directory_or_url, level=None):
+
+        self._directory_or_url = directory_or_url
+
+        self._is_url = is_url(directory_or_url)
+
+        self._properties = load_properties(directory_or_url)
+
+        self._tile_size = int(self._properties["hips_tile_width"])
+        self._order = int(self._properties["hips_order"])
+        self._level = self._order if level is None else level
+
+        self._tile_format = self._properties["hips_tile_format"]
+
+        self._nside = level_to_nside(self._level)
+        self._hp = HEALPix(nside=self._nside, frame="icrs", order="nested")
+        self._cdelt = 45 / self._tile_size / 2**self._level * 2**0.5
+
+        image_size = 5 * self._nside * self._tile_size
+
+        self.wcs = WCS(naxis=2)
+        self.wcs.wcs.ctype = "RA---HPX", "DEC--HPX"
+        self.wcs.wcs.crval = 0, 0
+        self.wcs.wcs.cdelt = self._cdelt, self._cdelt
+        self.wcs.wcs.crpix = image_size / 2, image_size / 2
+        self.wcs.wcs.crota = 0, 45
+        self.wcs.wcs.set()
+
+        self.shape = (image_size, image_size)
+        self.dtype = float
+        self.ndim = 2
+
+        self.chunksize = (self._tile_size, self._tile_size)
+
+        self._nan = np.nan * np.ones(self.chunksize, dtype=self.dtype)
+
+    def __getitem__(self, item):
+
+        # For now assume item is a list of slices. Find
+
+        imid = (item[0].start + item[0].stop) // 2
+        jmid = (item[1].start + item[1].stop) // 2
+
+        # Convert pixel coordinates to HEALPix indices
+
+        index = self._hp.skycoord_to_healpix(self.wcs.pixel_to_world(imid, jmid))
+
+        if index == -1:
+            return self._nan
+
+        return self._get_tile(level=self._level, index=index)
+
+    def _get_tile(self, *, level, index):
+
+        filename = tile_filename(
+            level=self._level,
+            index=index,
+            output_directory=self._directory_or_url,
+            extension="fits",
+        )
+
+        if self._is_url:
+            try:
+                return fits.getdata(filename).astype(float)[::-1]
+            except urllib.error.HTTPError:
+                return self._nan
+        else:
+            if os.path.exists(filename):
+                # FIXME: why flip vertically?
+                return fits.getdata(filename).astype(float)[::-1]
+            else:
+                return self._nan

reproject/hips/high_level.py

@@ -21,7 +21,9 @@
 from ..wcs_utils import has_celestial
 from .utils import (
     determine_healpix_level,
+    load_properties,
     make_tile_folders,
+    save_properties,
     tile_filename,
     tile_header,
 )
@@ -401,21 +403,6 @@ def save_index(directory):
         f.write(INDEX_HTML)
 
 
-def save_properties(directory, properties):
-    with open(os.path.join(directory, "properties"), "w") as f:
-        for key, value in properties.items():
-            f.write(f"{key:20s} = {value}\n")
-
-
-def load_properties(directory):
-    properties = {}
-    with open(os.path.join(directory, "properties")) as f:
-        for line in f:
-            key, value = line.split("=")
-            properties[key.strip()] = value.strip()
-    return properties
-
-
 def coadd_hips(input_directories, output_directory):
     """
     Given multiple HiPS directories, combine these into a single HiPS.

reproject/hips/utils.py

@@ -1,4 +1,5 @@
 import os
+import urllib
 
 import numpy as np
 from astropy.wcs import WCS
@@ -103,3 +104,31 @@ def determine_healpix_level(wcs_in, tile_size):
     target_level = nside_to_level(target_nside)
 
     return target_level
+
+
+def _is_url(directory):
+    return directory.startswith("http://") or directory.startswith("https://")
+
+
+def save_properties(directory, properties):
+    with open(os.path.join(directory, "properties"), "w") as f:
+        for key, value in properties.items():
+            f.write(f"{key:20s} = {value}\n")
+
+
+def load_properties(directory_or_url):
+
+    if _is_url(directory_or_url):
+        properties_filename, _ = urllib.request.urlretrieve(f"{directory_or_url}/properties")
+    else:
+        properties_filename = os.path.join(directory_or_url, "properties")
+
+    properties = {}
+    with open(properties_filename) as f:
+        for line in f:
+            if line.startswith("#"):
+                continue
+            key, value = line.split("=", 1)
+            properties[key.strip()] = value.strip()
+
+    return properties