Merge pull request #71 from imbasimba/min-max-cuts

Min, max, and suggested cuts

ci/azure-build-and-test.yml

@@ -5,36 +5,36 @@ parameters:
   type: object
   default:
 
-  - name: linux_37
-    vmImage: ubuntu-20.04
-    vars:
-      PYTHON_SERIES: "3.7"
-
   - name: linux_38
     vmImage: ubuntu-20.04
     vars:
       PYTHON_SERIES: "3.8"
 
-  - name: macos_37
-    vmImage: macos-10.15
+  - name: linux_39
+    vmImage: ubuntu-20.04
     vars:
-      PYTHON_SERIES: "3.7"
+      PYTHON_SERIES: "3.9"
 
   - name: macos_38
     vmImage: macos-10.15
     vars:
       PYTHON_SERIES: "3.8"
 
-  - name: windows_37
-    vmImage: windows-2019
+  - name: macos_39
+    vmImage: macos-10.15
     vars:
-      PYTHON_SERIES: "3.7"
+      PYTHON_SERIES: "3.9"
 
   - name: windows_38
     vmImage: windows-2019
     vars:
       PYTHON_SERIES: "3.8"
 
+  - name: windows_39
+    vmImage: windows-2019
+    vars:
+      PYTHON_SERIES: "3.9"
+
 jobs:
 - ${{ each build in parameters.builds }}:
   - job: ${{ format('build_{0}', build.name) }}

setup.py

@@ -78,7 +78,7 @@ def get_long_desc():
         "reproject",
         "shapely",
         "tqdm>=4.0",
-        "wwt_data_formats>=0.10.2",
+        "wwt_data_formats>=0.12",
     ],
     extras_require={
         "test": [

toasty/builder.py

@@ -190,7 +190,7 @@ def load_from_wwtl(self, cli_settings, wwtl_path, cli_progress=False):
         # existing imageset as much as possible, but update the parameters that
         # change in the tiling process.
 
-        wcs_keywords = self.imgset.wcs_headers_from_position()
+        wcs_keywords = self.imgset.wcs_headers_from_position(height=img.height)
         self.imgset.center_x = (
             self.imgset.center_y
         ) = 0  # hack to satisfy _check_no_wcs_yet()
@@ -225,6 +225,23 @@ def cascade(self, **kwargs):
         from .merge import averaging_merger, cascade_images
 
         cascade_images(self.pio, self.imgset.tile_levels, averaging_merger, **kwargs)
+        if "fits" in self.imgset.file_type:
+            from .pyramid import Pos
+            from astropy.io import fits
+            import numpy as np
+
+            with fits.open(
+                self.pio.tile_path(
+                    pos=Pos(n=0, x=0, y=0), format="fits", makedirs=False
+                )
+            ) as top_tile:
+                self.imgset.data_min = top_tile[0].header["DATAMIN"]
+                self.imgset.data_max = top_tile[0].header["DATAMAX"]
+                (
+                    self.imgset.pixel_cut_low,
+                    self.imgset.pixel_cut_high,
+                ) = np.nanpercentile(top_tile[0].data, [0.5, 99.5])
+
         return self
 
     def make_thumbnail_from_other(self, thumbnail_image):

toasty/fits_tiler.py

@@ -345,4 +345,11 @@ def _copy_hips_properties_to_builder(self):
         self.builder.imgset.base_degrees_per_tile = float(
             hips_properties["hips_initial_fov"]
         )
+        pixel_cut = hips_properties["hips_pixel_cut"].split(" ")
+        self.builder.imgset.pixel_cut_low = float(pixel_cut[0])
+        self.builder.imgset.pixel_cut_high = float(pixel_cut[1])
+        data_range = hips_properties["hips_data_range"].split(" ")
+        self.builder.imgset.data_min = float(data_range[0])
+        self.builder.imgset.data_max = float(data_range[1])
+
         self.builder.imgset.url = "Norder{0}/Dir{1}/Npix{2}"

toasty/image.py

@@ -634,8 +634,19 @@ def load_path(self, path):
                     wcs = WCS(hdul[0].header)
                 else:
                     wcs = None
-
-                img = Image.from_array(arr, wcs=wcs, default_format="fits")
+                max_value = self._get_header_value_or_none(
+                    header=hdul[0].header, keyword="DATAMAX"
+                )
+                min_value = self._get_header_value_or_none(
+                    header=hdul[0].header, keyword="DATAMIN"
+                )
+                img = Image.from_array(
+                    arr,
+                    wcs=wcs,
+                    default_format="fits",
+                    min_value=min_value,
+                    max_value=max_value,
+                )
             return img
 
         # Special handling for Photoshop files, used for some very large mosaics
@@ -680,6 +691,12 @@ def load_path(self, path):
         with open(path, "rb") as f:
             return self.load_stream(f)
 
+    def _get_header_value_or_none(self, header, keyword):
+        value = None
+        if keyword in header:
+            value = header[keyword]
+        return value
+
 
 class Image(object):
     """
@@ -694,6 +711,8 @@ class Image(object):
     _mode = None
     _default_format = "png"
     _wcs = None
+    _data_min = None
+    _data_max = None
 
     @classmethod
     def from_pil(cls, pil_img, wcs=None, default_format=None):
@@ -736,7 +755,9 @@ def from_pil(cls, pil_img, wcs=None, default_format=None):
         return inst
 
     @classmethod
-    def from_array(cls, array, wcs=None, default_format=None):
+    def from_array(
+        cls, array, wcs=None, default_format=None, min_value=None, max_value=None
+    ):
         """Create a new Image from an array-like data variable.
 
         Parameters
@@ -751,6 +772,12 @@ def from_array(cls, array, wcs=None, default_format=None):
             The default format to use when writing the image if none is
             specified explicitly. If not specified, this is automatically
             chosen at write time based on the array type.
+        min_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the lowest data value in this image and its children.
+        max_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the highest data value in this image and its children.
 
         Returns
         -------
@@ -778,6 +805,11 @@ def from_array(cls, array, wcs=None, default_format=None):
         inst._default_format = default_format or cls._default_format
         inst._array = array
         inst._wcs = wcs
+        if "fits" in inst._default_format:
+            if min_value is not None:
+                inst._data_min = min_value
+            if max_value is not None:
+                inst._data_max = max_value
         return inst
 
     def asarray(self):
@@ -871,6 +903,14 @@ def default_format(self, value):
         else:
             raise ValueError("Unrecognized format: {0}".format(value))
 
+    @property
+    def data_min(self):
+        return self._data_min
+
+    @property
+    def data_max(self):
+        return self._data_max
+
     def has_wcs(self):
         """
         Return whether this image has attached WCS information.
@@ -1075,7 +1115,9 @@ def update_into_maskable_buffer(self, buffer, iy_idx, ix_idx, by_idx, bx_idx):
                 f"unhandled mode `{self.mode}` in update_into_maskable_buffer"
             )
 
-    def save(self, path_or_stream, format=None, mode=None):
+    def save(
+        self, path_or_stream, format=None, mode=None, min_value=None, max_value=None
+    ):
         """
         Save this image to a filesystem path or stream
 
@@ -1084,6 +1126,18 @@ def save(self, path_or_stream, format=None, mode=None):
         path_or_stream : path-like object or file-like object
             The destination into which the data should be written. If file-like,
             the stream should accept bytes.
+        format : :class:`str` or ``None`` (the default)
+            The format name; one of ``SUPPORTED_FORMATS``
+        mode : :class:`toasty.image.ImageMode` or ``None`` (the default)
+            The image data mode to use if ``format`` is a ``PIL_FORMATS``
+        min_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the lowest data value in this image and its children.
+            If not set, the minimum value will be extracted from this image.
+        max_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the highest data value in this image and its children.
+            If not set, the maximum value will be extracted from this image.
         """
 
         _validate_format("format", format)
@@ -1107,14 +1161,20 @@ def save(self, path_or_stream, format=None, mode=None):
             # Avoid annoying RuntimeWarnings on all-NaN data
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
-
-                m = np.nanmin(arr)
-                if np.isfinite(m):  # Astropy will raise an error if we don't NaN-guard
-                    header["DATAMIN"] = m
-
-                m = np.nanmax(arr)
-                if np.isfinite(m):
-                    header["DATAMAX"] = m
+                if min_value is not None:
+                    header["DATAMIN"] = min_value
+                else:
+                    m = np.nanmin(arr)
+                    if np.isfinite(
+                        m
+                    ):  # Astropy will raise an error if we don't NaN-guard
+                        header["DATAMIN"] = m
+                if max_value is not None:
+                    header["DATAMAX"] = max_value
+                else:
+                    m = np.nanmax(arr)
+                    if np.isfinite(m):
+                        header["DATAMAX"] = m
 
             fits.writeto(
                 path_or_stream,

toasty/merge.py

@@ -165,13 +165,46 @@ def _cascade_images_serial(pio, start, merger, cli_progress):
                     )
 
             merged = Image.from_array(merger(buf.asarray()))
-            pio.write_image(pos, merged)
+            min_value, max_value = _get_min_max_of_children(
+                pio, [img0, img1, img2, img3]
+            )
+
+            pio.write_image(pos, merged, min_value=min_value, max_value=max_value)
             progress.update(1)
 
     if cli_progress:
         print()
 
 
+def _get_min_max_of_children(pio, children):
+    min_value = None
+    max_value = None
+    if "fits" in pio.get_default_format():
+        min_values = _get_existing_min_values(children)
+        if min_values:  # Check there are any valid min values
+            min_value = min(min_values)
+        max_values = _get_existing_max_values(children)
+        if max_values:  # Check there are any valid max values
+            max_value = max(max_values)
+    return min_value, max_value
+
+
+def _get_existing_min_values(images):
+    values = []
+    for image in images:
+        if image is not None and image.data_min is not None:
+            values.append(image.data_min)
+    return values
+
+
+def _get_existing_max_values(images):
+    values = []
+    for image in images:
+        if image is not None and image.data_max is not None:
+            values.append(image.data_max)
+    return values
+
+
 def _cascade_images_parallel(pio, start, merger, cli_progress, parallel):
     """Parallelized cascade operation
 
@@ -317,6 +350,9 @@ def _mp_cascade_worker(done_queue, ready_queue, done_event, pio, merger):
                     )
 
             merged = Image.from_array(merger(buf.asarray()))
-            pio.write_image(pos, merged)
+            min_value, max_value = _get_min_max_of_children(
+                pio, [img0, img1, img2, img3]
+            )
+            pio.write_image(pos, merged, min_value=min_value, max_value=max_value)
 
         done_queue.put(pos)

toasty/pyramid.py

@@ -338,7 +338,9 @@ def read_image(self, pos, default="none", masked_mode=None, format=None):
 
         return img
 
-    def write_image(self, pos, image, format=None, mode=None):
+    def write_image(
+        self, pos, image, format=None, mode=None, min_value=None, max_value=None
+    ):
         """Write an Image for the specified tile position.
 
         Parameters
@@ -347,10 +349,28 @@ def write_image(self, pos, image, format=None, mode=None):
             The tile position to write.
         image : :class:`toasty.image.Image`
             The image to write.
+        format : :class:`str` or ``None`` (the default)
+            The format name; one of ``SUPPORTED_FORMATS``
+        mode : :class:`toasty.image.ImageMode` or ``None`` (the default)
+            The image data mode to use if ``format`` is a ``PIL_FORMATS``
+        min_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the lowest data value in this image and its children.
+            If not set, the minimum value will be extracted from this image.
+        max_value : number or ``None`` (the default)
+            An optional number only used for FITS images.
+            The value represents to the highest data value in this image and its children.
+            If not set, the maximum value will be extracted from this image.
 
         """
         p = self.tile_path(pos, format=format or self._default_format)
-        image.save(p, format=format or self._default_format, mode=mode)
+        image.save(
+            p,
+            format=format or self._default_format,
+            mode=mode,
+            min_value=min_value,
+            max_value=max_value,
+        )
 
     @contextmanager
     def update_image(self, pos, default="none", masked_mode=None, format=None):