astropy · astrofrog · Oct 16, 2025 · Oct 15, 2025 · Oct 15, 2025 · Oct 15, 2025
diff --git a/CHANGES.md b/CHANGES.md
@@ -282,23 +282,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
 
@@ -307,34 +307,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
@@ -343,7 +343,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
 
@@ -352,20 +352,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
 
@@ -376,33 +376,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
 
@@ -413,22 +413,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
 

diff --git a/reproject/hips/_dask_array.py b/reproject/hips/_dask_array.py
@@ -12,6 +12,7 @@
 from astropy_healpix import HEALPix, level_to_nside
 from dask import array as da
 
+from ._trim_utils import fits_getdata_untrimmed
 from .high_level import VALID_COORD_SYSTEM
 from .utils import (
     is_url,
@@ -209,11 +210,14 @@ def _get_tile(self, *, level, index):
         else:
             filename = filename_or_url
 
-        with fits.open(filename) as hdulist:
-            hdu = hdulist[0]
-            data = hdu.data
-
-        return data
+        if self.ndim == 2:
+            return fits.getdata(filename)
+        else:
+            return fits_getdata_untrimmed(
+                filename,
+                tile_size=self._tile_width,
+                tile_depth=self._tile_depth,
+            )
 
 
 def hips_as_dask_array(directory_or_url, *, level=None):

diff --git a/reproject/hips/_trim_utils.py b/reproject/hips/_trim_utils.py
@@ -0,0 +1,51 @@
+import numpy as np
+from astropy.io import fits
+
+__all__ = ["fits_writeto_withtrim", "fits_getdata_untrimmed"]
+
+
+def fits_writeto_withtrim(filename, array, header, **kwargs):
+
+    if array.ndim == 3:
+
+        mask = ~np.isnan(array)
+
+        slices = []
+        for axis in range(3):
+            mask_1d = np.any(mask, axis=tuple(a for a in range(3) if a != axis))
+            indices = np.nonzero(mask_1d)[0]
+            if len(indices) == 0:
+                return
+            imin, imax = indices[0], indices[-1]
+            slices.append(slice(imin, imax + 1))
+            fits_axis = 3 - axis
+            header[f"TRIM{fits_axis}"] = imin
+            header[f"CRPIX{fits_axis}"] -= imin
+            header[f"NAXIS{fits_axis}"] = imax + 1 - imin
+        array = array[tuple(slices)]
+
+    array = array.astype(np.float32)
+    header["BITPIX"] = -32
+
+    fits.writeto(filename, array, header, **kwargs)
+
+
+def fits_getdata_untrimmed(filename, *, tile_size, tile_depth):
+
+    with fits.open(filename) as hdulist:
+        data = hdulist[0].data
+        header = hdulist[0].header
+
+    pad_before = tuple(header.get(f"TRIM{3 - axis}", 0) for axis in range(3))
+    shape = data.shape
+    pad_after = (
+        tile_depth - shape[0] - pad_before[0],
+        tile_size - shape[1] - pad_before[1],
+        tile_size - shape[2] - pad_before[2],
+    )
+
+    data = np.pad(data, list(zip(pad_before, pad_after, strict=False)))
+
+    assert data.shape == (tile_depth, tile_size, tile_size)
+
+    return data