Fixed a bug where obs coordinates were being interpreted as being on …

…cell centers instead of cell vertices. (#19)

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`.
 

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)

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",

pyproject.toml

@@ -16,7 +16,7 @@ classifiers = [
 dependencies = [
     "astropy",
     "requests",
-    "named-arrays==0.14.2",
+    "named-arrays==0.15.0",
 ]
 dynamic = ["version"]