First pass at a vbi_reproject function

[Cadair]

The objective of this PR in the short term is to make it easy to reproject a VBI mosaic to glue it together.

Here is a current example script:

import matplotlib.pyplot as plt
import reproject

import dkist
from dkist.processing.mosaicking import reproject_vbi  # New shiny thing

ds = dkist.load_dataset("~/dkist_data/pid_2_114/AGJZM/")

fig = plt.figure(figsize=(12, 12))

for i, tile in enumerate([d[0] for d in ds.flat]):
    ax = fig.add_subplot(ds.shape[0], ds.shape[1], i+1, projection=tile.wcs)
    ax.set_title(f"MINDEX1={tile.headers[0]['MINDEX1']}, MINDEX2={tile.headers[0]['MINDEX2']}")
    ax.imshow(tile.data)

fig.tight_layout()

new_ds, footprint = reproject_vbi(ds, reproject_function=reproject.reproject_interp, shape_out=(10900, 10900))

fig = plt.figure()
new_ds.plot()

[codecov]

Codecov Report

Attention: Patch coverage is 27.65957% with 34 lines in your changes are missing coverage. Please review.

Project coverage is 96.12%. Comparing base (788bfd9) to head (6512393).
Report is 1 commits behind head on main.

Files	Patch %	Lines
dkist/processing/mosaicking.py	27.65%	34 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #343      +/-   ##
==========================================
- Coverage   97.73%   96.12%   -1.62%     
==========================================
  Files          34       35       +1     
  Lines        1990     2037      +47     
==========================================
+ Hits         1945     1958      +13     
- Misses         45       79      +34     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[SolarDrew]

From a little basic testing so far, the helper function is not resilient against different kinds of WCS. For example, resampling the data before trying to reproject so that it can run in a sensible amount of time just breaks it:

In [6]: ds = dkist.load_dataset('/home/drew/sunpy/data/VBI_L1_20231016T184424_AGJZM.asdf')

In [7]: newtiles = []

In [8]: for i, tile in enumerate([d for d in ds.flat]):
   ...:     newtiles.append(tile.rebin((1, 4, 4), operation=np.sum))
   ...: 

In [9]: ds = dkist.TiledDataset(np.array(newtiles).reshape(ds.shape), inventory=ds.inventory)

In [10]: new_ds, footprint = reproject_vbi(ds, reproject_function=reproject.reproject_interp)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [10], line 1
----> 1 new_ds, footprint = reproject_vbi(ds, reproject_function=reproject.reproject_interp)

File ~/oss-projects/dkist/dkist/processing/mosaicking.py:64, in reproject_vbi(tds, edge_crop, reference_tile, reproject_function, roundtrip_coords, combine_function, shape_out)
     54 def reproject_vbi(
     55     tds: dkist.TiledDataset,
     56     *,
   (...)
     62     shape_out: tuple[int] = None,
     63 ):
---> 64     uni_tds = _unify_output_frames(tds, reference_tile)
     65     cropped = dkist.TiledDataset(
     66         np.array([ds[:, edge_crop:-edge_crop, edge_crop:-edge_crop] for ds in uni_tds.flat]).reshape(tds.shape),
     67         inventory=tds.inventory,
     68     )
     70     target_shape = shape_out or np.array(cropped[reference_tile][0].data.shape) * cropped.shape

File ~/oss-projects/dkist/dkist/processing/mosaicking.py:27, in _unify_output_frames(tds, reference_tile)
     22 """
     23 Given a `dkist.TiledDataset` return a new version where all the WCSes share
     24 the same celestial output frame.
     25 """
     26 ref_wcs = tds[reference_tile].wcs
---> 27 ref_celestial_frame = [f for f in ref_wcs.output_frame.frames if isinstance(f, cf.CelestialFrame)][0]
     28 new_datasets = []
     29 for ds in tds.flat:

AttributeError: 'HighLevelWCSWrapper' object has no attribute 'output_frame'