Fix wcs submap issues v2

tilemaker/metadata/fits.py

@@ -38,8 +38,17 @@
             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,17 +67,17 @@
                 )
 
             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]:
@@ -82,8 +91,8 @@
         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)
 

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,8 +140,8 @@
     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))
 
@@ -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

tilemaker/processing/test_submap_wcs.py

@@ -0,0 +1,112 @@
+# test_submap_wcs.py
+from astropy.wcs import WCS
+
+from tilemaker.processing.wcs_utils import build_submap_wcs
+
+# --- Inputs matching your real request ---
+LEFT = -72.3514
+RIGHT = -60.9358
+TOP = -39.4649
+BOTTOM = -41.4070
+
+TOLERANCE_DEG = 0.05
+
+
+def get_base_wcs() -> WCS:
+    """Reconstruct the same base WCS the layer provider returns."""
+    CDELT_RA = -0.0083333333333333
+    CDELT_DEC = 0.0083333333333333
+    NAXIS1 = int(360.0 / abs(CDELT_RA))
+    NAXIS2 = int(180.0 / abs(CDELT_DEC))
+
+    return WCS(
+        {
+            "NAXIS": 2,
+            "CRPIX1": NAXIS1 * 0.5,
+            "CRPIX2": NAXIS2 * 0.5 + 0.5,
+            "CRVAL1": 0.0,
+            "CRVAL2": 0.0,
+            "NAXIS1": NAXIS1,
+            "NAXIS2": NAXIS2,
+            "CDELT1": CDELT_RA,
+            "CDELT2": -CDELT_DEC,
+            "CTYPE1": "RA---CAR",
+            "CTYPE2": "DEC--CAR",
+            "CUNIT1": "deg",
+            "CUNIT2": "deg",
+            "LONPOLE": 0.0,
+            "LATPOLE": 90.0,
+            "RADESYS": "ICRS",
+        }
+    )
+
+
+def check(name, got, expected, tol=TOLERANCE_DEG):
+    diff = abs(got - expected)
+    status = "PASS" if diff < tol else "FAIL"
+    print(f"{status} {name}: got={got:.5f}  expected={expected:.5f}  diff={diff:.5f}")
+    return status == "PASS"
+
+
+def run():
+    base_wcs = get_base_wcs()
+    submap_wcs = build_submap_wcs(LEFT, RIGHT, TOP, BOTTOM, base_wcs)
+    padded_x_size, padded_y_size = submap_wcs.pixel_shape
+    x_size, y_size = submap_wcs.data_shape
+    offset_x = submap_wcs.data_offset_x
+    offset_y = submap_wcs.data_offset_y
+
+    print(f"\nPadded array size: x={padded_x_size}, y={padded_y_size}")
+    print(
+        f"Actual data size:  x={x_size}, y={y_size}  "
+        f"(offset_x={offset_x}, offset_y={offset_y})"
+    )
+    print(f"CRVAL:  {submap_wcs.wcs.crval}")
+    print(f"CRPIX:  {submap_wcs.wcs.crpix}")
+    print(f"CDELT:  {submap_wcs.wcs.cdelt}")
+    print()
+
+    # pixel_to_world uses 0-indexed pixels. Both axes are padded (see
+    # build_submap_wcs), so the actual cutout data lives at columns
+    # [offset_x, offset_x + x_size) and rows [offset_y, offset_y + y_size)
+    # rather than starting at pixel (0, 0).
+    #
+    # Which specific pixel corner maps to which specific world corner
+    # (e.g. does column offset_x hold LEFT or RIGHT?) is decided by
+    # base_wcs's own pixel-index-vs-world-value convention -- verified
+    # (via the real base layer's own tile-serving behavior, including a
+    # gradient marker to rule out mirroring) to need to match base_wcs's
+    # convention exactly, rather than a convention build_submap_wcs
+    # enforces itself. So instead of asserting a specific corner mapping,
+    # check the four corners as a set: two should be at RA=LEFT, two at
+    # RA=RIGHT (one each for the two rows), and two at Dec=TOP, two at
+    # Dec=BOTTOM (one each for the two columns) -- i.e. a correctly
+    # positioned, non-mirrored, non-rotated rectangle.
+    corners = [
+        submap_wcs.pixel_to_world(offset_x, offset_y),
+        submap_wcs.pixel_to_world(offset_x + x_size - 1, offset_y),
+        submap_wcs.pixel_to_world(offset_x, offset_y + y_size - 1),
+        submap_wcs.pixel_to_world(offset_x + x_size - 1, offset_y + y_size - 1),
+    ]
+    ras = sorted(c.ra.deg for c in corners)
+    decs = sorted(c.dec.deg for c in corners)
+
+    all_pass = all(
+        [
+            check("min RA  (two corners)", ras[0], LEFT),
+            check("min RA  (two corners)", ras[1], LEFT),
+            check("max RA  (two corners)", ras[2], RIGHT),
+            check("max RA  (two corners)", ras[3], RIGHT),
+            check("min Dec (two corners)", decs[0], BOTTOM),
+            check("min Dec (two corners)", decs[1], BOTTOM),
+            check("max Dec (two corners)", decs[2], TOP),
+            check("max Dec (two corners)", decs[3], TOP),
+        ]
+    )
+
+    print()
+    print("ALL PASS" if all_pass else "SOME CHECKS FAILED")
+
+
+if __name__ == "__main__":
+    run()

tilemaker/processing/wcs_utils.py

@@ -0,0 +1,153 @@
+import astropy.units as u
+import numpy as np
+from astropy.coordinates import SkyCoord
+from astropy.wcs import WCS
+
+
+class _CoordResult:
+    """Minimal ra/dec container. SkyCoord's constructor always re-wraps RA
+    to astropy's default [0, 360) range, discarding any custom wrap_angle,
+    so it can't be used to carry a -180..180 RA value through untouched."""
+
+    __slots__ = ("ra", "dec")
+
+    def __init__(self, ra, dec):
+        self.ra = ra
+        self.dec = dec
+
+
+class _ClientConventionWCS(WCS):
+    """A WCS whose pixel_to_world reports RA in -180..180 (client
+    convention) rather than astropy's default 0..360, since that's the
+    convention build_submap_wcs's left/right/CRVAL inputs are given in."""
+
+    def pixel_to_world(self, *pixel_arrays):
+        coord = super().pixel_to_world(*pixel_arrays)
+        return _CoordResult(coord.ra.wrap_at(180 * u.deg), coord.dec)
+
+
+def build_submap_wcs(
+    left: float,
+    right: float,
+    top: float,
+    bottom: float,
+    base_wcs: WCS,
+) -> WCS:
+    """
+    Build a WCS for a submap cutout.
+
+    Parameters
+    ----------
+    left, right : float
+        RA bounds in degrees (client convention, e.g. -180 to 180)
+    top, bottom : float
+        Dec bounds in degrees
+    base_wcs : WCS
+        The full-sky base WCS from the layer provider
+
+    Returns
+    -------
+    submap_wcs : WCS
+        `submap_wcs.pixel_shape` is the padded array shape the header
+        describes (see below). `submap_wcs.data_shape` is the actual
+        (x_size, y_size) of the requested cutout, and
+        `data_offset_x`/`data_offset_y` is where that data should be
+        written within an array sized `pixel_shape`.
+    """
+    # base_wcs may come straight from a FITS header with extra axes
+    # (frequency, Stokes, ...); reduce to the 2D celestial sub-WCS so
+    # crpix/crval/cdelt below always have exactly 2 elements.
+    base_wcs = base_wcs.celestial
+
+    cdelt_ra = base_wcs.wcs.cdelt[0]
+    cdelt_dec = base_wcs.wcs.cdelt[1]
+    cdelt1_mag = abs(cdelt_ra)
+    cdelt2_mag = abs(cdelt_dec)
+
+    x_size = int(round(abs(right - left) / cdelt1_mag))
+    y_size = int(round(abs(top - bottom) / cdelt2_mag))
+
+    bottom_left = SkyCoord(ra=left * u.deg, dec=bottom * u.deg)
+    top_right = SkyCoord(ra=right * u.deg, dec=top * u.deg)
+
+    # The header uses base_wcs's own CRPIX/CRVAL/LONPOLE/LATPOLE/CDELT-signs
+    # completely unmodified (just re-anchored to CRVAL=(0, 0) -- see below),
+    # not any re-derived values. Verified empirically against the real
+    # full-sky layer's own tile-serving behavior (including with a gradient
+    # marker to rule out mirroring): the tile-serving path
+    # (`providers/fits.py::extract_patch_from_fits`, and the tile-index
+    # remap in `server/layers.py`) reads straight from the FITS header's
+    # own WCS, and only produces correct, level-independent, seamless
+    # results when CRPIX/CDELT/LONPOLE/LATPOLE match a real base layer's
+    # own values exactly -- re-deriving *any* of them (a synthetic
+    # NAXIS/2-style CRPIX, a minimized/re-anchored CRPIX, a different
+    # LONPOLE/LATPOLE) broke positioning beyond the lowest zoom level in
+    # ways that were each individually hard to predict. This is the one
+    # configuration that's actually been confirmed correct, even though it
+    # costs more padding than a "minimal" CRPIX choice would.
+    ref_wcs = base_wcs.deepcopy()
+    ref_wcs.wcs.crval = [0.0, 0.0]
+
+    ref_crpix1, ref_crpix2 = ref_wcs.wcs.crpix
+    ref_cdelt1, ref_cdelt2 = ref_wcs.wcs.cdelt
+
+    px_bl, py_bl = ref_wcs.world_to_pixel(bottom_left)
+    px_tr, py_tr = ref_wcs.world_to_pixel(top_right)
+    px_bl, py_bl, px_tr, py_tr = float(px_bl), float(py_bl), float(px_tr), float(py_tr)
+
+    # Which world corner ends up at the smaller pixel index isn't something
+    # this function decides -- it falls out of base_wcs's own convention
+    # (its CDELT signs combined with its LONPOLE/LATPOLE), and is left
+    # alone rather than forced into a fixed orientation.
+    #
+    # The X (RA) axis is padded out to the width a genuine 360-degree-wide
+    # sky would have at base_wcs's own pixel scale -- not just "enough to
+    # reach the cutout (or CRPIX)". This costs more disk space per export
+    # (the array's width no longer scales down for small or nearby
+    # cutouts), but it's what makes CRPIX1 land at its own exact midpoint,
+    # the same way it already does on base_wcs's own full-sky grid. That
+    # midpoint property is what the server's tile-index "flip" fold
+    # (server/layers.py::get_tile) and its per-tile mirror
+    # (processing/renderer.py) both assume for the RA axis; without it,
+    # "flip" scrambles a submap layer's tiles instead of leaving them
+    # alone. With it, a submap-derived layer needs no special-casing at
+    # all -- it's handled by the exact same code path as a directly
+    # registered full-sky FITS file.
+    #
+    # The Y (Dec) axis does NOT get the same full-height treatment: both
+    # the fold and the mirror only ever act on the RA axis (a Dec value
+    # doesn't have a "0-360 vs -180-180" convention to reconcile), so Y
+    # keeps the original minimal padding -- just enough to reach the
+    # cutout's own far edge or CRPIX2, whichever is further. Padding Y out
+    # to a full 180 degrees actively breaks things here: base_wcs's own
+    # CRPIX2 reflects wherever its real (often Dec-limited) survey data
+    # sits, not the midpoint of a true pole-to-pole span, so a
+    # full-height array pushes far pixel rows outside the CAR
+    # projection's valid range and get_bbox() ends up with NaN corners.
+    naxis1_padded = int(round(360.0 / cdelt1_mag))
+    naxis2_padded = int(np.ceil(max(py_bl, py_tr, ref_crpix2))) + 1
+    data_offset_x = int(round(min(px_bl, px_tr)))
+    data_offset_y = int(round(min(py_bl, py_tr)))
+
+    submap_wcs = _ClientConventionWCS(naxis=2)
+    submap_wcs.wcs.ctype = base_wcs.wcs.ctype
+    submap_wcs.wcs.cunit = base_wcs.wcs.cunit
+    submap_wcs.wcs.radesys = base_wcs.wcs.radesys
+    submap_wcs.wcs.lonpole = ref_wcs.wcs.lonpole
+    submap_wcs.wcs.latpole = ref_wcs.wcs.latpole
+
+    submap_wcs.wcs.crval = ref_wcs.wcs.crval
+    submap_wcs.wcs.crpix = [ref_crpix1, ref_crpix2]
+    submap_wcs.wcs.cdelt = [ref_cdelt1, ref_cdelt2]
+
+    # Carry the cutout size on the WCS itself (astropy's (NAXIS1, NAXIS2)
+    # pixel_shape convention) so callers can recover it without a second
+    # return value. Both axes are padded (see above); data_offset_x/y is
+    # where the actual x_size x y_size cutout data should be written within
+    # that wider/taller array.
+    submap_wcs.pixel_shape = (naxis1_padded, naxis2_padded)
+    submap_wcs.data_shape = (x_size, y_size)
+    submap_wcs.data_offset_x = data_offset_x
+    submap_wcs.data_offset_y = data_offset_y
+
+    return submap_wcs

tilemaker/providers/fits.py

@@ -261,7 +261,18 @@ def extract_patch_from_fits(
         log = log.bind(dt=end - start)
         log.debug("fits.no_data")
 
-        return None
+        # The requested window falls entirely outside the FITS array (as
+        # opposed to overlapping it but landing on NaN padding, which
+        # extract_array handles fine and doesn't hit this branch). Return an
+        # all-NaN patch of the same (post-subsample) shape a real cutout
+        # would have, rather than None: this keeps every "no data here"
+        # tile behaving the same way (a normal, cacheable blank tile
+        # response) regardless of which of the two cases produced it.
+        # Returning None instead surfaces as an HTTP 404 for this specific
+        # case only, which client tile-loading code can treat very
+        # differently from a successful-but-blank tile.
+        blank_shape = tuple(s // subsample_every for s in shape)
+        return np.full(blank_shape, np.nan)
 
     if subsample_every > 1:
         log = log.bind(subsample_every=subsample_every)