Fixes for HiPS functionality

CHANGES.md

@@ -296,23 +296,23 @@
 ## 0.9 - 2022-09-02
 
 - Drop support for Python 3.7.
-- 
+-
 - Infrastructure and packaging updates.
-- 
+-
 - Made many improvements, bug fixes, and significant speed-ups for the adaptive
 - resampling algorithm, `reproject_adaptive`. These bug fixes may cause
 - changes to the reprojected images, which are typically negligible.
 - Improvements include the addition of a flux-conserving mode, support for a
 - Gaussian filter kernel, a menu of boundary-handling modes, and a
 - `center_jacobian` flag to trade speed for accuracy with rapidly-varying
 - transformations.
-- 
+-
 - Added a `roundtrip_coords` argument to `reproject_adaptive` and
 - `reproject_interp`. By default, all coordinate transformations are run in
 - both directions to handle some situations where they are ambiguous. This can
 - be disabled by setting `roundtrip_coords=False` which may offer a
 - significant speed increase.
-- 
+-
 
 ## 0.8 - 2021-08-11
 
@@ -321,34 +321,34 @@
 - possible to solve for the optimal WCS for a set of images that couldn't fit
 - into memory all at once, since the actual data aren't needed for optimal WCS
 - determination. [#242]
-- 
+-
 - Fix implementation of `hdu_weights` in `reproject_and_coadd`. [#249]
-- 
+-
 
 ## 0.7.1 - 2020-05-29
 
 - Fixed compatibility with Astropy 4.1. [#234]
-- 
+-
 - Updated minimum requirement for SciPy. [#236]
-- 
+-
 
 ## 0.7 - 2020-04-02
 
 - Made C extension in overlapArea.c thread-safe by removing global
 - variables. [#211]
-- 
+-
 - Made it possible to control whether to output debugging information
 - from overlapArea.c by setting DEBUG_OVERLAP_AREA=1 at build-time. [#211]
-- 
+-
 - Fix compatibility with astropy v4.0.1. [#227]
-- 
+-
 - Disable parallelization by default in `reproject_exact` - this can be
 - enabled with `parallel=True`. [#227]
-- 
+-
 - Fixed a bug with `reproject_exact` with `parallel=False` and
 - `return_footprint=False`, which caused the footprint to be returned
 - anyway. [#227]
-- 
+-
 - The infrastructure of the package has been updated in line with the
 - APE 17 roadmap (https://github.com/astropy/astropy-APEs/blob/main/APE17.rst).
 - The main changes are that the `python setup.py test` and
@@ -357,7 +357,7 @@
 - (https://tox.readthedocs.io) package and run `tox -e test` and
 - `tox -e build_docs`. It is also possible to run pytest and sphinx
 - directly. [#228]
-- 
+-
 
 ## 0.6 - 2019-11-01
 
@@ -366,20 +366,20 @@
 - `independent_celestial_slices=` argument to `reproject_interp` has
 - been deprecated since it is no longer needed, as transformations are
 - automatically done in the most efficient way possible. [#166]
-- 
+-
 - Include a warning for high resolution images with `reproject_exact`,
 - since if the pixels are <0.05", precision issues can occur. [#200]
-- 
+-
 - Added a new `reproject_and_coadd` function for doing mosaicking of
 - individual images, and added section in documentation about mosaicking.
 - [#186]
-- 
+-
 - Added a new reproject.adaptive sub-package that implements the DeForest
 - (2004) algorithm for reprojection. [#52]
-- 
+-
 - Fixed a bug that caused 'exact' reprojection results to have numerical
 - issues when doing identity transformations. [#190]
-- 
+-
 
 ## 0.5.1 - 2019-09-01
 
@@ -390,33 +390,33 @@
 
 - Improve parse_output_projection to make it so that the output projection
 - can be specified as a filename. [#150]
-- 
+-
 - Fixed a bug that caused HEALPix maps in RING order to not be correctly
 - interpreted. [#163]
-- 
+-
 - Make it possible to specify the output array for reprojection using the
 - `output_array=` keyword argument. [#115]
-- 
+-
 
 ## 0.4 - 2018-01-29
 
 - Refactored HEALPix reprojection code to use the astropy-healpix package
 - instead of healpy. [#139]
-- 
+-
 - Added the ability to specify an output array in `reproject_interp`, which
 - permits the use of memory-mapped arrays and therefore provides the capability
 - to handle data cubes much larger than memory [#115]
-- 
+-
 - Fix test 32-bit test failures. [#146]
-- 
+-
 - Fix an issue with reprojecting images where there are two solutions along
 - the line of sight by forcing round-tripping of coordinate conversions [#129]
-- 
+-
 - Explicitly define default HDU as 0 for normal reprojection and 1 for
 - HEALPix reprojection. [#119]
-- 
+-
 - Added a function to find the optimal WCS for a set of images. [#136, #137]
-- 
+-
 
 ## 0.3.2 - 2017-10-22
 
@@ -427,22 +427,22 @@
 ## 0.3.1 - 2016-07-07
 
 - Include missing license file in tarball.
-- 
+-
 - Updated documentation to remove warnings about early versions.
-- 
+-
 
 ## 0.3 - 2016-07-06
 
 - Allow users to pass a `field=` option to `reproject_from_healpix`
 - to access different fields in a HEALPIX file. [#86]
-- 
+-
 - Significant improvements to performance when the input data is a large
 - memory-mapped array. [#105]
-- 
+-
 - Significant refactoring of interpolating reprojection to improve support for
 - n-dimensional arrays, optionally including two celestial axes (in which
 - case the coordinate transformation is taken into account). [#96, #102]
-- 
+-
 
 ## 0.2 - 2015-10-29
 

reproject/hips/_dask_array.py

@@ -30,11 +30,13 @@ class HiPSArray:
 
     def __init__(self, directory_or_url, level=None, level_depth=None):
 
-        self._directory_or_url = directory_or_url
+        # We strip any trailing slashes since we then assume in the rest of the
+        # code that we need to add a slash (and double slashes cause issues for URLs)
+        self._directory_or_url = str(directory_or_url).rstrip("/")
 
-        self._is_url = is_url(directory_or_url)
+        self._is_url = is_url(self._directory_or_url)
 
-        self._properties = load_properties(directory_or_url)
+        self._properties = load_properties(self._directory_or_url)
 
         if self._properties["dataproduct_type"] == "image":
             self.ndim = 2
@@ -51,7 +53,7 @@ def __init__(self, directory_or_url, level=None, level_depth=None):
         else:
             if level > self._order_spatial:
                 raise ValueError(
-                    f"HiPS dataset at {directory_or_url} does not contain spatial level {level} data"
+                    f"HiPS dataset at {self._directory_or_url} does not contain spatial level {level} data"
                 )
             elif level < 0:
                 raise ValueError("level should be positive")
@@ -70,7 +72,7 @@ def __init__(self, directory_or_url, level=None, level_depth=None):
             else:
                 if level_depth > self._order_depth:
                     raise ValueError(
-                        f"HiPS dataset at {directory_or_url} does not contain spectral level {level_depth} data"
+                        f"HiPS dataset at {self._directory_or_url} does not contain spectral level {level_depth} data"
                     )
                 elif level_depth < 0:
                     raise ValueError("level_depth should be positive")

reproject/hips/high_level.py

@@ -1,6 +1,7 @@
 import os
 import shutil
 import uuid
+import warnings
 from concurrent.futures import ThreadPoolExecutor
 from copy import deepcopy
 from datetime import datetime
@@ -470,7 +471,7 @@ def process(index):
     if tile_format == "fits":
         generated_properties["hips_pixel_bitpix"] = -64
         if not np.isinf(pixel_min) and not np.isinf(pixel_max):
-            properties["hips_pixel_cut"] = (pixel_min, pixel_max)
+            properties["hips_pixel_cut"] = f"{pixel_min} {pixel_max}"
 
     generated_properties.update(properties)
 
@@ -482,7 +483,6 @@ def process(index):
         output_directory=output_directory,
         ndim=ndim,
         frame=frame,
-        tile_dims=tile_dims,
         tile_format=tile_format,
         tile_size=tile_size,
         tile_depth=tile_depth,
@@ -522,7 +522,6 @@ def compute_lower_resolution_tiles(
     output_directory,
     ndim,
     frame,
-    tile_dims,
     tile_format,
     tile_size,
     tile_depth,
@@ -531,6 +530,11 @@ def compute_lower_resolution_tiles(
     indices=None,
 ):
 
+    if ndim == 3:
+        tile_dims = tile_size, tile_depth
+    else:
+        tile_dims = tile_size
+
     # Iterate over higher levels and compute lower resolution tiles
 
     half_tile_size = tile_size // 2
@@ -599,7 +603,10 @@ def compute_lower_resolution_tiles(
 
                         if tile_format == "fits":
                             tile_data = fits.getdata(subtile_filename)
-                            data = block_reduce(tile_data, 2, func=np.mean)
+                            # np.nanmean can emit warnings, so ignore these
+                            with warnings.catch_warnings():
+                                warnings.simplefilter("ignore")
+                                data = block_reduce(tile_data, 2, func=np.nanmean)
                         else:
                             data = block_reduce(
                                 np.array(Image.open(subtile_filename))[::-1],
@@ -633,13 +640,16 @@ def compute_lower_resolution_tiles(
 
                         if os.path.exists(subtile_filename):
 
-                            data = block_reduce(
-                                fits_getdata_untrimmed(
-                                    subtile_filename, tile_size=tile_size, tile_depth=tile_depth
-                                ),
-                                2,
-                                func=np.mean,
-                            )
+                            # np.nanmean can emit warnings, so ignore these
+                            with warnings.catch_warnings():
+                                warnings.simplefilter("ignore")
+                                data = block_reduce(
+                                    fits_getdata_untrimmed(
+                                        subtile_filename, tile_size=tile_size, tile_depth=tile_depth
+                                    ),
+                                    2,
+                                    func=np.nanmean,
+                                )
 
                             if subindex_spec == 0:
                                 subtile_slice = [slice(None, half_tile_depth)]

reproject/hips/utils.py

@@ -50,14 +50,14 @@ def spectral_index_to_coord(level, index):
     )
 
 
-def spectral_coord_to_index(level, coord):
-    return (
-        np.floor(
-            2 ** (level + 1)
-            * np.log10(coord.to_value(u.Hz) / FREQ_MIN)
-            / np.log10(FREQ_MAX / FREQ_MIN)
-        )
-    ).astype(int)
+def spectral_coord_to_index(level, coord, rounded=True):
+    index = (
+        2 ** (level + 1) * np.log10(coord.to_value(u.Hz) / FREQ_MIN) / np.log10(FREQ_MAX / FREQ_MIN)
+    )
+    if rounded:
+        return np.floor(index).astype(int)
+    else:
+        return index
 
 
 def map_header(*, level, frame, tile_dims):