From c1f9c5b55788d4d2d1486af6eb8dbc64785c446f Mon Sep 17 00:00:00 2001
From: Jeremy Myers <inernestmusic@gmail.com>
Date: Thu, 21 May 2026 12:47:24 -0400
Subject: [PATCH] Include WCS header in submap

---
 tilemaker/processing/extractor.py | 30 ++++++++++++++++++++++++++----
 tilemaker/server/layers.py        |  5 +++--
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/tilemaker/processing/extractor.py b/tilemaker/processing/extractor.py
index 50e476d..3a47a08 100644
--- a/tilemaker/processing/extractor.py
+++ b/tilemaker/processing/extractor.py
@@ -24,7 +24,7 @@ def extract(
     metadata: DataConfiguration,
     grants: set[str],
     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).
 
@@ -84,11 +84,11 @@ def extract(
     wcs = WCS(
         {
             "NAXIS": 2,
+            "NAXIS1": NAXIS1,
+            "NAXIS2": NAXIS2,
             "CRPIX1": NAXIS1 * 0.5,
             "CRPIX2": NAXIS2 * 0.5 + 0.5,
             "CRVAL1": 0.0,
-            "NAXIS1": NAXIS1,
-            "NAXIS2": NAXIS2,
             "CRVAL2": 0.0,
             "CDELT1": CDELT_RA,
             "CDELT2": -CDELT_DEC,
@@ -197,4 +197,26 @@ def extract(
 
     log = log.info("extractor.complete")
 
-    return buffer, pushables
+    # Create WCS for submap; note that only CRPIX and NAXIS change
+    submap_wcs = WCS(
+        {
+            "NAXIS": 2,
+            "NAXIS1": x_size,
+            "NAXIS2": y_size,
+            "CRPIX1": wcs.wcs.crpix[0] - left_pix,
+            "CRPIX2": wcs.wcs.crpix[1] - bottom_pix,
+            "CRVAL1": wcs.wcs.crval[0],
+            "CRVAL2": wcs.wcs.crval[1],
+            "CDELT1": wcs.wcs.cdelt[0],
+            "CDELT2": wcs.wcs.cdelt[1],
+            "CTYPE1": "RA---CAR",
+            "CTYPE2": "DEC--CAR",
+            "CUNIT1": "deg",
+            "CUNIT2": "deg",
+            "LONPOLE": 90.0,
+            "LATPOLE": 0.0,
+            "RADESYS": "ICRS",
+        }
+    )
+
+    return buffer, pushables, submap_wcs
diff --git a/tilemaker/server/layers.py b/tilemaker/server/layers.py
index d6dbead..7f72df2 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, wcs = extract(
         layer_id=layer_id,
         left=left,
         right=right,
@@ -186,7 +186,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 = wcs.to_header()
+            hdu = fits.PrimaryHDU(submap, header=header)
             hdu.writeto(output)
             return Response(content=output.getvalue(), media_type="image/fits")