From a49490277eb15e05ca4776a2afbd31f350b60f6d Mon Sep 17 00:00:00 2001 From: Jeremy Myers Date: Mon, 27 Jul 2026 12:57:56 -0400 Subject: [PATCH 1/3] Fix WCS issues with extracted submaps --- tilemaker/metadata/fits.py | 29 +++-- tilemaker/processing/extractor.py | 29 ++++- tilemaker/processing/test_submap_wcs.py | 112 ++++++++++++++++++++ tilemaker/processing/wcs_utils.py | 135 ++++++++++++++++++++++++ tilemaker/providers/fits.py | 13 ++- tilemaker/server/layers.py | 39 ++++++- 6 files changed, 340 insertions(+), 17 deletions(-) create mode 100644 tilemaker/processing/test_submap_wcs.py create mode 100644 tilemaker/processing/wcs_utils.py diff --git a/tilemaker/metadata/fits.py b/tilemaker/metadata/fits.py index 9a5e025..ad35377 100644 --- a/tilemaker/metadata/fits.py +++ b/tilemaker/metadata/fits.py @@ -38,8 +38,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,17 +67,17 @@ 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]: 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..15365e4 --- /dev/null +++ b/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() diff --git a/tilemaker/processing/wcs_utils.py b/tilemaker/processing/wcs_utils.py new file mode 100644 index 0000000..c9c689e --- /dev/null +++ b/tilemaker/processing/wcs_utils.py @@ -0,0 +1,135 @@ +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 written array starts at pixel 0 of base_wcs's own grid and spans + # out to whichever pixel index is furthest from it -- either the + # cutout's own far edge, or (if base_wcs's CRPIX itself sits beyond + # the cutout) that reference point, since `extract_patch_from_fits` + # treats `abs(CRPIX1) > NAXIS1` as a double-wrap sign and "corrects" + # it, corrupting the projection if the array doesn't reach that far. + naxis1_padded = int(np.ceil(max(px_bl, px_tr, ref_crpix1))) + 1 + 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 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..9f00429 100644 --- a/tilemaker/server/layers.py +++ b/tilemaker/server/layers.py @@ -158,7 +158,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 +168,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) @@ -186,7 +187,8 @@ def get_submap( return Response(content=output.getvalue(), media_type="image/png") elif ext == "fits": with io.BytesIO() as output: - hdu = fits.PrimaryHDU(submap) + header = submap_wcs.to_header() + hdu = fits.PrimaryHDU(submap, header) hdu.writeto(output) return Response(content=output.getvalue(), media_type="image/fits") @@ -242,8 +244,24 @@ 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. - if level != 0: + # it's a card-folding operation. This is only meaningful for a layer + # whose own pixel grid spans the full sky (RA=0 sits at the exact + # horizontal midpoint of its array) -- true for a directly-registered + # full-sky FITS file, but not for a submap cutout (see + # processing/extractor.py / processing/wcs_utils.py), whose array is + # much narrower and whose CRPIX1 does not sit at its midpoint. + # Applying the fold to such a layer scrambles tile positions instead + # of leaving them alone, so it's gated on the layer's own bounding + # box actually spanning (close to) 360 degrees of RA. + layer = next( + (lyr for lyr in request.app.config.layers if lyr.layer_id == layer_id), + None, + ) + spans_full_sky = layer is not None and ( + abs(layer.bounding_right - layer.bounding_left) > 350 + ) + + if spans_full_sky and level != 0: # Level of zero requires no flipping apart from at the tile level. midpoint = 2 ** (level) if x < midpoint: @@ -251,6 +269,19 @@ def get_tile( else: x = (2 ** (level) - 1) - (x - midpoint) + midpoint + if not spans_full_sky: + # renderer.render()'s own per-tile mirror is the other half of + # this same full-sky-only "flip" mechanism (it's what the above + # index remap is meant to be paired with). Left enabled here, + # it mirrors each tile in isolation around its own center -- + # for a full-sky layer that's fine since every tile's content + # is part of one continuous card-folded whole, but for a + # submap cutout it pushes each tile's real data away from the + # tile it's adjacent to, opening a visible gap wherever real + # data straddles a tile boundary. Since this layer doesn't + # need the fold at all, just don't mirror its tiles. + render_options.flip = False + if ext not in ["jpg", "webp", "png"]: raise HTTPException(status_code=400, detail="Not an acceptable extension") From d4673f51b0e164dc8491697d7dd676e947e8ca2f Mon Sep 17 00:00:00 2001 From: Jeremy Myers Date: Tue, 28 Jul 2026 10:45:54 -0400 Subject: [PATCH 2/3] Fix WCS submap issues --- tilemaker/processing/wcs_utils.py | 32 ++++++++++---- tilemaker/server/layers.py | 69 +++++++++++++++---------------- 2 files changed, 58 insertions(+), 43 deletions(-) diff --git a/tilemaker/processing/wcs_utils.py b/tilemaker/processing/wcs_utils.py index c9c689e..e070f68 100644 --- a/tilemaker/processing/wcs_utils.py +++ b/tilemaker/processing/wcs_utils.py @@ -100,13 +100,31 @@ def build_submap_wcs( # (its CDELT signs combined with its LONPOLE/LATPOLE), and is left # alone rather than forced into a fixed orientation. # - # The written array starts at pixel 0 of base_wcs's own grid and spans - # out to whichever pixel index is furthest from it -- either the - # cutout's own far edge, or (if base_wcs's CRPIX itself sits beyond - # the cutout) that reference point, since `extract_patch_from_fits` - # treats `abs(CRPIX1) > NAXIS1` as a double-wrap sign and "corrects" - # it, corrupting the projection if the array doesn't reach that far. - naxis1_padded = int(np.ceil(max(px_bl, px_tr, ref_crpix1))) + 1 + # 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))) diff --git a/tilemaker/server/layers.py b/tilemaker/server/layers.py index 9f00429..2362812 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, @@ -186,11 +189,28 @@ 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: - header = submap_wcs.to_header() - hdu = fits.PrimaryHDU(submap, header) - hdu.writeto(output) - return Response(content=output.getvalue(), media_type="image/fits") + # A submap-derived layer's array is now padded out to a full-sky + # -sized grid (see processing/wcs_utils.py::build_submap_wcs), so + # `submap` can be multiple GB even though almost all of it is NaN + # padding. Serializing through an in-memory io.BytesIO -- and then + # Response(content=...) taking another copy via output.getvalue() + # -- means holding several multiples of that size in memory at + # once, which can exhaust memory outright for a large base layer. + # Writing directly to a file and streaming it back via + # FileResponse instead keeps astropy's write side to its own + # internal (small, chunked) buffering, and avoids the extra + # in-memory copies entirely. + 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( @@ -244,24 +264,14 @@ 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. This is only meaningful for a layer - # whose own pixel grid spans the full sky (RA=0 sits at the exact - # horizontal midpoint of its array) -- true for a directly-registered - # full-sky FITS file, but not for a submap cutout (see - # processing/extractor.py / processing/wcs_utils.py), whose array is - # much narrower and whose CRPIX1 does not sit at its midpoint. - # Applying the fold to such a layer scrambles tile positions instead - # of leaving them alone, so it's gated on the layer's own bounding - # box actually spanning (close to) 360 degrees of RA. - layer = next( - (lyr for lyr in request.app.config.layers if lyr.layer_id == layer_id), - None, - ) - spans_full_sky = layer is not None and ( - abs(layer.bounding_right - layer.bounding_left) > 350 - ) - - if spans_full_sky and level != 0: + # it's a card-folding operation. This assumes the layer's own pixel + # grid spans the full sky with RA=0 at the exact horizontal + # midpoint of its array -- true for a directly-registered full-sky + # FITS file, and also true for a submap cutout (see + # processing/wcs_utils.py::build_submap_wcs), whose array is + # padded out to a full-sky-sized grid specifically so this holds + # for it too, rather than needing to be special-cased here. + if level != 0: # Level of zero requires no flipping apart from at the tile level. midpoint = 2 ** (level) if x < midpoint: @@ -269,19 +279,6 @@ def get_tile( else: x = (2 ** (level) - 1) - (x - midpoint) + midpoint - if not spans_full_sky: - # renderer.render()'s own per-tile mirror is the other half of - # this same full-sky-only "flip" mechanism (it's what the above - # index remap is meant to be paired with). Left enabled here, - # it mirrors each tile in isolation around its own center -- - # for a full-sky layer that's fine since every tile's content - # is part of one continuous card-folded whole, but for a - # submap cutout it pushes each tile's real data away from the - # tile it's adjacent to, opening a visible gap wherever real - # data straddles a tile boundary. Since this layer doesn't - # need the fold at all, just don't mirror its tiles. - render_options.flip = False - if ext not in ["jpg", "webp", "png"]: raise HTTPException(status_code=400, detail="Not an acceptable extension") From 815648ae62df746d2f75aba0c29c94fdf5c09544 Mon Sep 17 00:00:00 2001 From: Jeremy Myers Date: Wed, 12 Aug 2026 12:54:20 -0400 Subject: [PATCH 3/3] Refactor WCS fix to reduce file size --- pyproject.toml | 3 + tilemaker/metadata/fits.py | 65 ++--- tilemaker/processing/test_submap_wcs.py | 307 ++++++++++++++++++------ tilemaker/processing/wcs_utils.py | 131 ++++++---- tilemaker/server/layers.py | 34 ++- 5 files changed, 374 insertions(+), 166 deletions(-) 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 ad35377..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.""" @@ -82,38 +117,10 @@ def sanitize_nonscalar(x): 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/test_submap_wcs.py b/tilemaker/processing/test_submap_wcs.py index 15365e4..ece42e8 100644 --- a/tilemaker/processing/test_submap_wcs.py +++ b/tilemaker/processing/test_submap_wcs.py @@ -1,24 +1,37 @@ -# test_submap_wcs.py -from astropy.wcs import WCS +""" +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 -from tilemaker.processing.wcs_utils import build_submap_wcs +import astropy.units as u +import numpy as np +import pytest +import structlog +from astropy.io import fits +from astropy.wcs import WCS -# --- Inputs matching your real request --- -LEFT = -72.3514 -RIGHT = -60.9358 -TOP = -39.4649 -BOTTOM = -41.4070 +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 -TOLERANCE_DEG = 0.05 +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 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)) +def _base_wcs() -> WCS: return WCS( { "NAXIS": 2, @@ -28,8 +41,8 @@ def get_base_wcs() -> WCS: "CRVAL2": 0.0, "NAXIS1": NAXIS1, "NAXIS2": NAXIS2, - "CDELT1": CDELT_RA, - "CDELT2": -CDELT_DEC, + "CDELT1": -CDELT, + "CDELT2": CDELT, "CTYPE1": "RA---CAR", "CTYPE2": "DEC--CAR", "CUNIT1": "deg", @@ -41,47 +54,24 @@ def get_base_wcs() -> WCS: ) -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" +# 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 run(): - base_wcs = get_base_wcs() +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) - 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 + offset_x, offset_y = submap_wcs.data_offset_x, 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), @@ -91,22 +81,203 @@ def run(): 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), - ] + # 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] ) - print() - print("ALL PASS" if all_pass else "SOME CHECKS FAILED") + 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 + ) -if __name__ == "__main__": - run() + 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 index e070f68..1fb89e1 100644 --- a/tilemaker/processing/wcs_utils.py +++ b/tilemaker/processing/wcs_utils.py @@ -1,8 +1,18 @@ +import math + import astropy.units as u -import numpy as np 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 @@ -64,27 +74,27 @@ def build_submap_wcs( 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)) + 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. 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. + # 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] @@ -99,35 +109,51 @@ def build_submap_wcs( # 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. # - # 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))) + # 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 @@ -137,17 +163,20 @@ def build_submap_wcs( 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.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 (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) + # 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 = data_offset_x - submap_wcs.data_offset_y = data_offset_y + 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/server/layers.py b/tilemaker/server/layers.py index 2362812..7e3f7e7 100644 --- a/tilemaker/server/layers.py +++ b/tilemaker/server/layers.py @@ -189,17 +189,15 @@ def get_submap( renderer.render(output, submap, render_options=render_options) return Response(content=output.getvalue(), media_type="image/png") elif ext == "fits": - # A submap-derived layer's array is now padded out to a full-sky - # -sized grid (see processing/wcs_utils.py::build_submap_wcs), so - # `submap` can be multiple GB even though almost all of it is NaN - # padding. Serializing through an in-memory io.BytesIO -- and then - # Response(content=...) taking another copy via output.getvalue() - # -- means holding several multiples of that size in memory at - # once, which can exhaust memory outright for a large base layer. - # Writing directly to a file and streaming it back via - # FileResponse instead keeps astropy's write side to its own - # internal (small, chunked) buffering, and avoids the extra - # in-memory copies entirely. + # 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) @@ -264,13 +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. This assumes the layer's own pixel - # grid spans the full sky with RA=0 at the exact horizontal - # midpoint of its array -- true for a directly-registered full-sky - # FITS file, and also true for a submap cutout (see - # processing/wcs_utils.py::build_submap_wcs), whose array is - # padded out to a full-sky-sized grid specifically so this holds - # for it too, rather than needing to be special-cased here. + # 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)