diff --git a/pyproject.toml b/pyproject.toml index f508e39..6453c9b 100644 --- a/pyproject.toml +++ b/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" diff --git a/tilemaker/metadata/fits.py b/tilemaker/metadata/fits.py index 9a5e025..b0c066e 100644 --- a/tilemaker/metadata/fits.py +++ b/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.""" diff --git a/tilemaker/processing/extractor.py b/tilemaker/processing/extractor.py index 50e476d..85720fb 100644 --- a/tilemaker/processing/extractor.py +++ b/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 @@ def extract( 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 @@ def extract( 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 @@ def extract( } ) + 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") @@ -197,4 +206,20 @@ def extract( 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 diff --git a/tilemaker/processing/test_submap_wcs.py b/tilemaker/processing/test_submap_wcs.py new file mode 100644 index 0000000..ece42e8 --- /dev/null +++ b/tilemaker/processing/test_submap_wcs.py @@ -0,0 +1,283 @@ +""" +Tests for build_submap_wcs, covering both correct positioning of the +cutout and the stride-aligned padding scheme that replaced padding the RA +axis out to a full-sky-sized array (see processing/wcs_utils.py). +""" + +import tempfile +from pathlib import Path + +import astropy.units as u +import numpy as np +import pytest +import structlog +from astropy.io import fits +from astropy.wcs import WCS + +from tilemaker.metadata.fits import FITSLayerProvider, tile_size_for_scale +from tilemaker.processing.wcs_utils import _PADDING_SAFETY_MARGIN, build_submap_wcs +from tilemaker.providers.core import PullableTile +from tilemaker.providers.fits import FITSTileProvider, extract_patch_from_fits + +LOG = structlog.get_logger() + +# A base layer at ~0.176 deg/pixel: NAXIS1=2048, NAXIS2=1024, which lands on +# the "clean 256" branch of tile_size_for_scale and yields a 4-level +# pyramid (levels 0-3, coarsest subsample stride = 2**(4-1) = 8) -- enough +# levels to meaningfully exercise stride-phase alignment while staying +# small/fast. +NAXIS1 = 2048 +NAXIS2 = 1024 +CDELT = 360.0 / NAXIS1 + + +def _base_wcs() -> WCS: + 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, + "CDELT2": CDELT, + "CTYPE1": "RA---CAR", + "CTYPE2": "DEC--CAR", + "CUNIT1": "deg", + "CUNIT2": "deg", + "LONPOLE": 0.0, + "LATPOLE": 90.0, + "RADESYS": "ICRS", + } + ) + + +# A small cutout close to (RA, Dec) = (180, -90), which for this base WCS +# (CRVAL=(0,0), CRPIX at the array midpoint) lands close to native pixel +# (0, 0) -- keeps the synthetic "base layer" array small in the tile-serving +# regression test below while still exercising a real, non-trivial pixel +# offset (not exactly zero). +LEFT, RIGHT = 176.0, 179.0 +BOTTOM, TOP = -89.0, -86.0 + + +def test_build_submap_wcs_places_corners_correctly(): + """The four corners of the cutout should land at (LEFT/RIGHT, TOP/BOTTOM), + not mirrored or rotated.""" + base_wcs = _base_wcs() + submap_wcs = build_submap_wcs(LEFT, RIGHT, TOP, BOTTOM, base_wcs) + + x_size, y_size = submap_wcs.data_shape + offset_x, offset_y = submap_wcs.data_offset_x, submap_wcs.data_offset_y + + 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) + + # left/right/top/bottom get rounded to the nearest whole pixel (see + # x_size/y_size in build_submap_wcs), so allow slack of a couple of + # pixels at this test's coarse (~0.18 deg/pixel) synthetic resolution. + tol = 2 * CDELT + assert ras[0] == pytest.approx(LEFT, abs=tol) + assert ras[1] == pytest.approx(LEFT, abs=tol) + assert ras[2] == pytest.approx(RIGHT, abs=tol) + assert ras[3] == pytest.approx(RIGHT, abs=tol) + assert decs[0] == pytest.approx(BOTTOM, abs=tol) + assert decs[1] == pytest.approx(BOTTOM, abs=tol) + assert decs[2] == pytest.approx(TOP, abs=tol) + assert decs[3] == pytest.approx(TOP, abs=tol) + + +def test_build_submap_wcs_padding_is_stride_aligned_not_full_sky(): + """The whole point of the fix: padding should be a small, bounded + multiple of the tile pyramid's coarsest subsample stride, not a + full-360-degree-wide array.""" + base_wcs = _base_wcs() + submap_wcs = build_submap_wcs(LEFT, RIGHT, TOP, BOTTOM, base_wcs) + + _, number_of_levels = tile_size_for_scale(CDELT * u.deg, CDELT * u.deg) + stride = 2 ** (number_of_levels - 1) + pad = _PADDING_SAFETY_MARGIN * stride + + naxis1, naxis2 = submap_wcs.pixel_shape + x_size, y_size = submap_wcs.data_shape + + # Padded array must be much smaller than a full-sky-sized array would + # have been (the old behavior padded NAXIS1 to 360 / cdelt ~= NAXIS1). + assert naxis1 < NAXIS1 // 4 + assert naxis2 < NAXIS2 // 4 + + # But still a whole multiple of `pad` (so CRPIX lands on a stride- + # aligned boundary), and big enough to hold the actual cutout. + assert naxis1 % pad == 0 + assert naxis2 % pad == 0 + assert naxis1 >= x_size + assert naxis2 >= y_size + + # CRPIX must have shifted by a whole multiple of `pad` relative to the + # base layer's own CRPIX -- that's what keeps subsample phase aligned + # with the source layer's own grid at every zoom level. + base_crpix1, base_crpix2 = base_wcs.wcs.crpix + shift_x = base_crpix1 - submap_wcs.wcs.crpix[0] + shift_y = base_crpix2 - submap_wcs.wcs.crpix[1] + assert shift_x % pad == 0 + assert shift_y % pad == 0 + + +def _flip_tile_x(x: int, level: int) -> int: + """Mirrors the tile-index fold in server/layers.py::get_tile.""" + if level == 0: + return x + midpoint = 2**level + if x < midpoint: + return (2**level - 1) - x + return (2**level - 1) - (x - midpoint) + midpoint + + +def test_submap_tile_serving_matches_base_layer_pixels(): + """Re-ingesting a submap export and serving tiles from it (at any zoom + level, with or without `flip`) must reproduce the same native pixel + values the *source* layer would have served for the same sky location. + + This is a regression test for the stride-phase bug found during + investigation: a naive "minimal CRPIX" cutout reproduces the correct + pixels only at the finest zoom level, and silently reads the wrong + native pixels (a phase-shifted subsample) at every coarser level. + """ + base_wcs = _base_wcs() + submap_wcs = build_submap_wcs(LEFT, RIGHT, TOP, BOTTOM, base_wcs) + + naxis1, naxis2 = submap_wcs.pixel_shape + x_size, y_size = submap_wcs.data_shape + data_offset_x, data_offset_y = submap_wcs.data_offset_x, submap_wcs.data_offset_y + + base_crpix1, base_crpix2 = base_wcs.wcs.crpix + aligned_offset_x = round(base_crpix1 - submap_wcs.wcs.crpix[0]) + aligned_offset_y = round(base_crpix2 - submap_wcs.wcs.crpix[1]) + + # A lightweight stand-in for the base layer's real array: only the + # small local footprint the padded submap actually touches (this test's + # LEFT/RIGHT/BOTTOM/TOP were chosen close to native pixel (0, 0) so this + # stays small), encoded with each pixel's own native (x, y) index so any + # phase shift in subsampling shows up as a wrong value rather than + # coincidentally matching. + margin = 4 + footprint_x = aligned_offset_x + naxis1 + margin + footprint_y = aligned_offset_y + naxis2 + margin + assert footprint_x < 500 and footprint_y < 500, ( + "test cutout constants no longer land near native pixel (0, 0); " + "the synthetic base array would be unexpectedly large" + ) + + base_data = ( + np.arange(footprint_x, dtype=np.float64)[None, :] * 100_000 + + np.arange(footprint_y, dtype=np.float64)[:, None] + ) + + submap_data = np.full((naxis2, naxis1), np.nan) + submap_data[ + data_offset_y : data_offset_y + y_size, data_offset_x : data_offset_x + x_size + ] = base_data[ + aligned_offset_y + data_offset_y : aligned_offset_y + data_offset_y + y_size, + aligned_offset_x + data_offset_x : aligned_offset_x + data_offset_x + x_size, + ] + + base_header = base_wcs.to_header() + submap_header = submap_wcs.to_header() + + _, number_of_levels = tile_size_for_scale(CDELT * u.deg, CDELT * u.deg) + provider = FITSTileProvider.__new__(FITSTileProvider) + + mismatches = [] + tested_any = False + for flip in (False, True): + for level in range(number_of_levels): + subsample_every = 2 ** (number_of_levels - level - 1) + n_tiles_x = 2 ** (level + 1) + n_tiles_y = 2**level + for display_x in range(n_tiles_x): + for y in range(n_tiles_y): + fetch_x = _flip_tile_x(display_x, level) if flip else display_x + info = provider._get_tile_info( + PullableTile( + layer_id="test", x=fetch_x, y=y, level=level, grants=None + ) + ) + + # Only test tiles that overlap the cutout at all -- + # elsewhere both sides are legitimately all-NaN/absent. + ra0, ra1 = sorted(info["ra_range"]) + dec0, dec1 = sorted(info["dec_range"]) + if ra1 < LEFT or ra0 > RIGHT or dec1 < BOTTOM or dec0 > TOP: + continue + + base_hdu = fits.PrimaryHDU(base_data, header=base_header.copy()) + submap_hdu = fits.PrimaryHDU( + submap_data, header=submap_header.copy() + ) + + base_patch = extract_patch_from_fits( + hdu=base_hdu, subsample_every=subsample_every, log=LOG, **info + ) + submap_patch = extract_patch_from_fits( + hdu=submap_hdu, + subsample_every=subsample_every, + log=LOG, + **info, + ) + tested_any = True + + if base_patch.shape != submap_patch.shape: + mismatches.append((flip, level, display_x, y, "shape mismatch")) + continue + + both_finite = np.isfinite(base_patch) & np.isfinite(submap_patch) + if not both_finite.any(): + continue + if not np.array_equal( + base_patch[both_finite], submap_patch[both_finite] + ): + mismatches.append((flip, level, display_x, y, "value mismatch")) + + assert tested_any, "no tiles overlapped the test cutout -- test setup is broken" + assert not mismatches, f"{len(mismatches)} tile mismatch(es): {mismatches[:10]}" + + +def test_get_bbox_recovers_tight_bounds_not_whole_sky(): + """A re-ingested submap layer's reported bounding box should be close + to the requested cutout, not the entire sky -- a bug the old full-sky + padding caused as a side effect (CRPIX left unmodified inside an array + padded out to a full-sky width made the corner pixels of the *padded* + array read out as +/-180 RA, +/-90 Dec).""" + base_wcs = _base_wcs() + submap_wcs = build_submap_wcs(LEFT, RIGHT, TOP, BOTTOM, base_wcs) + + naxis1, naxis2 = submap_wcs.pixel_shape + data = np.zeros((naxis2, naxis1), dtype=np.float64) + header = submap_wcs.to_header() + hdu = fits.PrimaryHDU(data, header=header) + + with tempfile.TemporaryDirectory() as tmp: + path = Path(tmp) / "submap.fits" + hdu.writeto(path) + + provider = FITSLayerProvider(filename=path) + bbox = provider.get_bbox() + + pad_deg = ( + _PADDING_SAFETY_MARGIN + * 2 ** (tile_size_for_scale(CDELT * u.deg, CDELT * u.deg)[1] - 1) + * CDELT + ) + + assert bbox["bounding_left"] == pytest.approx(LEFT, abs=pad_deg + 0.1) + assert bbox["bounding_right"] == pytest.approx(RIGHT, abs=pad_deg + 0.1) + assert bbox["bounding_top"] == pytest.approx(TOP, abs=pad_deg + 0.1) + assert bbox["bounding_bottom"] == pytest.approx(BOTTOM, abs=pad_deg + 0.1) diff --git a/tilemaker/processing/wcs_utils.py b/tilemaker/processing/wcs_utils.py new file mode 100644 index 0000000..1fb89e1 --- /dev/null +++ b/tilemaker/processing/wcs_utils.py @@ -0,0 +1,182 @@ +import math + +import astropy.units as u +from astropy.coordinates import SkyCoord +from astropy.wcs import WCS + +from tilemaker.metadata.fits import tile_size_for_scale + +# How many multiples of the tile pyramid's coarsest subsample stride to pad +# by, beyond the bare minimum needed to keep subsample phase aligned (see +# the comment in build_submap_wcs below). A larger margin costs a bit more +# file size but absorbs edge effects at the very coarsest zoom levels for +# large/off-center cutouts. +_PADDING_SAFETY_MARGIN = 4 + + +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 = round(abs(right - left) / cdelt1_mag) + y_size = 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, other than a whole-pixel shift applied + # below to both CRPIX and the array origin together. 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`) + # 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 modulo a + # whole multiple of the tile pyramid's coarsest subsample stride (see + # below) -- shifting by anything else broke positioning beyond the + # lowest zoom level in ways that were each individually hard to + # predict, because `extract_patch_from_fits`'s subsampling picks its + # phase from the *absolute* pixel position of each request. + 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. + raw_offset_x = round(min(px_bl, px_tr)) + raw_offset_y = round(min(py_bl, py_tr)) + + # A submap FITS file gets re-ingested as a brand-new layer via + # `tilemaker open`, whose own tile pyramid depth is re-derived from + # scratch purely from pixel scale (FITSLayerProvider.calculate_tile_size, + # which shares tile_size_for_scale with this function) -- so it always + # matches what we compute here, independent of the *source* layer's own + # (possibly config-overridden) number_of_levels. + # + # FITSTileProvider.pull() requests native pixels at subsample stride + # 2**(number_of_levels - 1 - level), whose coarsest value is + # 2**(number_of_levels - 1). extract_array/overlap_slices pick which + # native pixels land in a given subsampled tile based on the *absolute* + # pixel position of the request, so shifting CRPIX (and the array's + # pixel-0 origin) by anything that isn't a whole multiple of that + # coarsest stride changes which native pixels a coarse-zoom tile reads + # -- scrambling positioning at every level except the finest. Padding + # to a whole multiple of the stride (with a small safety margin, + # _PADDING_SAFETY_MARGIN) keeps the same subsample phase as base_wcs's + # own (unshifted) grid at every level, without needing to pad all the + # way out to a full-sky-sized array. This applies identically to both + # axes -- the tile-index "flip" fold (server/layers.py::get_tile) and + # its per-tile mirror (processing/renderer.py) only remap an abstract + # tile index derived from a hardcoded virtual full-sky grid + # (FITSTileProvider._get_tile_info); they never read CRPIX/NAXIS, so + # they need no special-casing here at all. + _, number_of_levels = tile_size_for_scale(cdelt1_mag * u.deg, cdelt2_mag * u.deg) + stride = max(1, 2 ** (number_of_levels - 1)) + pad = _PADDING_SAFETY_MARGIN * stride + + # Defensive upper bound: never pad past the size a genuine full-sky + # array would have at this pixel scale (what today's export used + # unconditionally) -- this keeps aligned_far_* sane even if `pad` were + # ever unexpectedly large relative to the cutout. + max_naxis_x = round(360.0 / cdelt1_mag) + max_naxis_y = round(180.0 / cdelt2_mag) + + aligned_offset_x = max(0, (raw_offset_x // pad) * pad) + aligned_offset_y = max(0, (raw_offset_y // pad) * pad) + aligned_far_x = min(max_naxis_x, math.ceil((raw_offset_x + x_size) / pad) * pad) + aligned_far_y = min(max_naxis_y, math.ceil((raw_offset_y + y_size) / pad) * pad) + + naxis1 = aligned_far_x - aligned_offset_x + naxis2 = aligned_far_y - aligned_offset_y + + 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 - aligned_offset_x, + ref_crpix2 - aligned_offset_y, + ] + 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 out to a stride-aligned boundary + # (see above); data_offset_x/y is where the actual x_size x y_size + # cutout data should be written within that array. + submap_wcs.pixel_shape = (naxis1, naxis2) + submap_wcs.data_shape = (x_size, y_size) + submap_wcs.data_offset_x = raw_offset_x - aligned_offset_x + submap_wcs.data_offset_y = raw_offset_y - aligned_offset_y + + return submap_wcs diff --git a/tilemaker/providers/fits.py b/tilemaker/providers/fits.py index d976eb4..8b03e34 100644 --- a/tilemaker/providers/fits.py +++ b/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) diff --git a/tilemaker/server/layers.py b/tilemaker/server/layers.py index d6dbead..7e3f7e7 100644 --- a/tilemaker/server/layers.py +++ b/tilemaker/server/layers.py @@ -3,6 +3,8 @@ """ import io +import os +import tempfile from typing import Literal from astropy.io import fits @@ -14,6 +16,7 @@ Request, Response, ) +from fastapi.responses import FileResponse from tilemaker.metadata.definitions import ( BandMenuState, @@ -158,7 +161,7 @@ def get_submap( Get a submap of the specified band. """ - submap, pushables = extract( + submap, pushables, submap_wcs = extract( layer_id=layer_id, left=left, right=right, @@ -168,6 +171,7 @@ def get_submap( grants=request.auth.scopes, metadata=request.app.config, show_grid=show_grid, + is_fits=ext == "fits", ) bt.add_task(request.app.tiles.push, pushables) @@ -185,10 +189,26 @@ def get_submap( renderer.render(output, submap, render_options=render_options) return Response(content=output.getvalue(), media_type="image/png") elif ext == "fits": - with io.BytesIO() as output: - hdu = fits.PrimaryHDU(submap) - hdu.writeto(output) - return Response(content=output.getvalue(), media_type="image/fits") + # submap's array is padded slightly beyond the requested cutout + # (see processing/wcs_utils.py::build_submap_wcs) to keep a + # re-ingested copy's tile-serving subsample phase aligned with the + # source layer's own -- a small, bounded amount, not a full-sky + # -sized grid. Still, write directly to a file and stream it back + # via FileResponse rather than through an in-memory io.BytesIO (and + # the extra copy Response(content=...) would take via + # output.getvalue()), so astropy's write side keeps to its own + # internal (small, chunked) buffering. + header = submap_wcs.to_header() + hdu = fits.PrimaryHDU(submap, header) + tmp = tempfile.NamedTemporaryFile(suffix=".fits", delete=False) + tmp.close() + hdu.writeto(tmp.name, overwrite=True) + bt.add_task(os.remove, tmp.name) + return FileResponse( + tmp.name, + media_type="image/fits", + filename=f"{layer_id}_submap.fits", + ) def core_tile_retrieval( @@ -242,7 +262,13 @@ def get_tile( if render_options.flip: # Flipping is really a reconfiguration of -180 < RA < 180 to 360 < RA < 0; - # it's a card-folding operation. + # it's a card-folding operation on the abstract tile index (x, + # level), derived from a hardcoded virtual full-sky grid + # (providers/fits.py::FITSTileProvider._get_tile_info). It never + # reads the underlying file's own WCS/CRPIX/NAXIS, so it works + # identically for a directly-registered full-sky FITS file and for + # a re-ingested submap cutout (processing/wcs_utils.py:: + # build_submap_wcs) without any special-casing here. if level != 0: # Level of zero requires no flipping apart from at the tile level. midpoint = 2 ** (level)