Fix wcs submap issues v2

pyproject.toml

@@ -30,6 +30,9 @@ dependencies = [
     "astropydantic"
 ]
 
+[project.optional-dependencies]
+dev = ["pytest"]
+
 [project.scripts]
 tilemaker = "tilemaker.client.cli:main"
 tilemaker-db = "tilemaker.metadata.database:main"

tilemaker/metadata/fits.py

@@ -14,6 +14,41 @@
 from pydantic import BaseModel
 
 
+def tile_size_for_scale(scale_x_deg, scale_y_deg) -> tuple[int, int]:
+    """
+    Derive (tile_size, number_of_levels) for a tile pyramid covering the
+    full sky at the given per-axis pixel scale. Shared between
+    FITSLayerProvider.calculate_tile_size (for a directly-ingested file) and
+    processing.wcs_utils.build_submap_wcs (which needs to know, ahead of
+    ingestion, what number_of_levels a from-scratch `tilemaker open` of its
+    export would derive) so the two always agree.
+    """
+    # The full sky spans 360 deg in RA, 180 deg in Dec
+    map_size_x = math.floor(360 * units.deg / scale_x_deg)
+    map_size_y = math.floor(180 * units.deg / scale_y_deg)
+
+    max_size = max(map_size_x, map_size_y)
+
+    # See if 256 fits.
+    if (map_size_x % 256 == 0) and (map_size_y % 256 == 0):
+        tile_size = 256
+        number_of_levels = int(math.log2(max_size // 256))
+        return tile_size, number_of_levels
+
+    # Oh no, remove all the powers of two until
+    # we get an odd number.
+    this_tile_size = map_size_y
+
+    # Also don't make it too small.
+    while this_tile_size % 2 == 0 and this_tile_size > 512:
+        this_tile_size = this_tile_size // 2
+
+    number_of_levels = int(math.log2(max_size // this_tile_size))
+    tile_size = this_tile_size
+
+    return tile_size, number_of_levels
+
+
 class LayerProvider(BaseModel):
     """Base class for layer providers."""
 
@@ -38,8 +73,17 @@ def get_bbox(self) -> dict[str, float]:
             data = handle[self.hdu]
             wcs = WCS(header=data.header)
 
-            top_right = wcs.array_index_to_world(*[0] * data.header.get("NAXIS", 2))
-            bottom_left = wcs.array_index_to_world(*[x - 1 for x in data.data.shape])
+            # Evaluate both opposite corners and take min/max explicitly,
+            # rather than assuming pixel (0, 0) is always the max-RA/max-Dec
+            # corner: that assumption holds for a typical telescope FITS
+            # file (RA decreasing, Dec increasing with pixel index) but not
+            # in general -- e.g. a submap cut out via
+            # tilemaker.processing.wcs_utils.build_submap_wcs can have
+            # either axis running the other way, since its orientation is
+            # derived from how the pixel buffer was actually assembled, not
+            # from this convention.
+            corner_a = wcs.array_index_to_world(*[0] * data.header.get("NAXIS", 2))
+            corner_b = wcs.array_index_to_world(*[x - 1 for x in data.data.shape])
 
             def sanitize(x):
                 return (
@@ -58,53 +102,25 @@ def sanitize_nonscalar(x):
                 )
 
             try:
-                tr = sanitize(top_right)
-                bl = sanitize(bottom_left)
+                a = sanitize(corner_a)
+                b = sanitize(corner_b)
             except TypeError:
-                tr = sanitize_nonscalar(top_right)
-                bl = sanitize_nonscalar(bottom_left)
+                a = sanitize_nonscalar(corner_a)
+                b = sanitize_nonscalar(corner_b)
 
         return {
-            "bounding_left": bl[0].value,
-            "bounding_right": tr[0].value,
-            "bounding_top": tr[1].value,
-            "bounding_bottom": bl[1].value,
+            "bounding_left": min(a[0], b[0]).value,
+            "bounding_right": max(a[0], b[0]).value,
+            "bounding_top": max(a[1], b[1]).value,
+            "bounding_bottom": min(a[1], b[1]).value,
         }
 
     def calculate_tile_size(self) -> tuple[int, int]:
         """Calculate appropriate tile size based on FITS file properties."""
-        # Need to figure out how big the whole 'map' is, i.e. moving it up
-        # so that it fills the whole space.
         wcs = self.get_wcs()
 
         scale = wcs.proj_plane_pixel_scales()
-        scale_x_deg = scale[0]
-        scale_y_deg = scale[1]
-
-        # The full sky spans 360 deg in RA, 180 deg in Dec
-        map_size_x = int(math.floor(360 * units.deg / scale_x_deg))
-        map_size_y = int(math.floor(180 * units.deg / scale_y_deg))
-
-        max_size = max(map_size_x, map_size_y)
-
-        # See if 256 fits.
-        if (map_size_x % 256 == 0) and (map_size_y % 256 == 0):
-            tile_size = 256
-            number_of_levels = int(math.log2(max_size // 256))
-            return tile_size, number_of_levels
-
-        # Oh no, remove all the powers of two until
-        # we get an odd number.
-        this_tile_size = map_size_y
-
-        # Also don't make it too small.
-        while this_tile_size % 2 == 0 and this_tile_size > 512:
-            this_tile_size = this_tile_size // 2
-
-        number_of_levels = int(math.log2(max_size // this_tile_size))
-        tile_size = this_tile_size
-
-        return tile_size, number_of_levels
+        return tile_size_for_scale(scale[0], scale[1])
 
     def get_wcs(self) -> WCS:
         """Get the WCS object from the FITS file."""

tilemaker/processing/extractor.py

@@ -11,6 +11,7 @@
 from astropy.wcs import WCS
 
 from tilemaker.metadata.core import DataConfiguration
+from tilemaker.processing.wcs_utils import build_submap_wcs
 from tilemaker.providers.core import PullableTile, PushableTile, Tiles
 
 
@@ -23,8 +24,9 @@
     tiles: Tiles,
     metadata: DataConfiguration,
     grants: set[str],
+    is_fits: bool,
     show_grid: bool = False,
-) -> tuple[np.array, list[PushableTile]]:
+) -> tuple[np.array, list[PushableTile], WCS]:
     """
     Extract a sub-map from a band between RA and Dec ranges (in degrees).
 
@@ -46,6 +48,8 @@
         Metadata object
     grants: set[str]
         Grants of the requesting user
+    is_fits: bool
+        Used to determine whether or not we need to derive a submap_wcs
     show_grid: bool = False
         Whether to 'show' the grid (grids are set as NaN values)
     """
@@ -102,6 +106,11 @@
         }
     )
 
+    submap_wcs = None
+
+    if is_fits:
+        submap_wcs = build_submap_wcs(left, right, top, bottom, base_wcs)
+
     # Convert RA/Dec to pixel values. No idea why we need to take the negative here.
     # Probably something I don't understand about wcs.
     tr = SkyCoord(ra=right, dec=top, unit="deg")
@@ -131,10 +140,10 @@
     buffer = np.zeros((int(y_size), (int(x_size))))
 
     # Figure out which tiles we overlap.
     end_tile_x = int(math.ceil(float(right_pix) / layer.tile_size))
     start_tile_x = int(math.floor(float(left_pix) / layer.tile_size))
     end_tile_y = int(math.ceil(top_pix / layer.tile_size))
     start_tile_y = int(math.floor(bottom_pix / layer.tile_size))
 
     # Load the tiles and push the data into the buffer.
     pushables = []
@@ -197,4 +206,20 @@
 
     log = log.info("extractor.complete")
 
-    return buffer, pushables
+    if is_fits:
+        # submap_wcs's axes are padded (see build_submap_wcs) so that its
+        # CRPIX stays within NAXIS on both axes -- write the tight buffer
+        # into a larger, NaN-filled array at the offsets submap_wcs
+        # actually describes, rather than writing it out at its own tight
+        # size. This only affects the FITS export; PNG/JPG/WEBP renders
+        # still use the tight buffer directly, unaffected.
+        padded_x_size, padded_y_size = submap_wcs.pixel_shape
+        offset_x = submap_wcs.data_offset_x
+        offset_y = submap_wcs.data_offset_y
+        padded_buffer = np.full((padded_y_size, padded_x_size), np.nan)
+        padded_buffer[
+            offset_y : offset_y + int(y_size), offset_x : offset_x + int(x_size)
+        ] = buffer
+        buffer = padded_buffer
+
+    return buffer, pushables, submap_wcs