From c1fc2f0be582dbea124cbc2c90007e779e331518 Mon Sep 17 00:00:00 2001
From: Roy Smart <roytsmart@gmail.com>
Date: Mon, 4 Nov 2024 13:25:04 -0700
Subject: [PATCH] Fixed a bug where obs coordinates were being interpreted as
 being on cell centers instead of cell vertices. (#19)

---
 docs/index.rst                    |  2 +-
 iris/sg/_spectrograph.py          |  2 +-
 iris/sg/background/_background.py | 21 ++++++++++++++++-----
 pyproject.toml                    |  2 +-
 4 files changed, 19 insertions(+), 8 deletions(-)

diff --git a/docs/index.rst b/docs/index.rst
index d835355..dd7ae88 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -1,7 +1,7 @@
 Introduction
 ============
 
-The `Interface Region Imaging Spectograph <iris.lmsal.com>`_ (IRIS) is a NASA
+The `Interface Region Imaging Spectograph <https://iris.lmsal.com>`_ (IRIS) is a NASA
 Small Explorer satellite which has been taking continuous ultraviolet images of
 the Sun since 2013 :cite:p:`DePontieu2014`.
 
diff --git a/iris/sg/_spectrograph.py b/iris/sg/_spectrograph.py
index 8f0ad72..41abe07 100644
--- a/iris/sg/_spectrograph.py
+++ b/iris/sg/_spectrograph.py
@@ -399,7 +399,7 @@ def empty(
                     ),
                 ),
             ),
-            shape_wcs=shape_wcs,
+            shape_wcs={a: shape_wcs[a] + 1 for a in shape_wcs},
         )
 
         shape = na.broadcast_shapes(shape_base, shape_wcs)
diff --git a/iris/sg/background/_background.py b/iris/sg/background/_background.py
index 21c80e3..cf2c6d2 100644
--- a/iris/sg/background/_background.py
+++ b/iris/sg/background/_background.py
@@ -73,7 +73,7 @@ def average(
             )
     """
     obs = obs.copy_shallow()
-    shape = obs.shape
+    shape = obs.inputs.shape
     obs.inputs = na.TemporalSpectralPositionalVectorArray(
         time=obs.inputs.time.ndarray.jd.mean(),
         wavelength=obs.inputs.wavelength.broadcast_to(shape).mean(axis),
@@ -313,8 +313,15 @@ def fit(
         # brightest spectral line
         wavelength_center = avg.wavelength_center.ndarray.mean()
 
+        # Interpolate wavelength samples onto cell centers
+        wavelength = avg.inputs.wavelength
+        for a in wavelength.shape:
+            lower = {a: slice(None, ~0)}
+            upper = {a: slice(+1, None)}
+            wavelength = (wavelength[lower] + wavelength[upper]) / 2
+
         # Convert wavelength to velocity units
-        velocity = avg.inputs.wavelength.to(
+        velocity = wavelength.to(
             unit=u.km / u.s,
             equivalencies=u.doppler_optical(wavelength_center),
         )
@@ -507,7 +514,11 @@ def subtract_spectral_line(
 
     where_crop = np.isfinite(obs.outputs).mean(obs.axis_wavelength) > 0.7
 
-    velocity = obs.inputs.wavelength[where_crop]
+    velocity = obs.inputs.wavelength
+    for a in velocity.shape:
+        velocity = (velocity[{a: slice(None, ~0)}] + velocity[{a: slice(1, None)}]) / 2
+
+    velocity = velocity[where_crop]
 
     where = np.abs(velocity) < 150 * u.km / u.s
 
@@ -781,12 +792,12 @@ def estimate(
         # Plot the result
         with astropy.visualization.quantity_support():
             fig, ax = plt.subplots()
-            na.plt.plot(
+            na.plt.stairs(
                 obs.inputs.wavelength.mean(obs.axis_time),
                 np.nanmedian(obs.outputs, axis=axis_txy),
                 label="original",
             )
-            na.plt.plot(
+            na.plt.stairs(
                 obs_nobg.inputs.wavelength.mean(obs.axis_time),
                 np.nanmedian(obs_nobg.outputs, axis=axis_txy),
                 label="corrected",
diff --git a/pyproject.toml b/pyproject.toml
index 26ef34a..3e2eab5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -16,7 +16,7 @@ classifiers = [
 dependencies = [
     "astropy",
     "requests",
-    "named-arrays==0.14.2",
+    "named-arrays==0.15.0",
 ]
 dynamic = ["version"]