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..e070f68 --- /dev/null +++ b/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 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..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, @@ -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,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: - hdu = fits.PrimaryHDU(submap) - 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( @@ -242,7 +264,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. 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)