Added Offset example and tweaked times

docs/guide.rst

@@ -122,13 +122,10 @@ We use the following command to read and load the data from a SJI level 2 file:
     Observatory:           IRIS
     Instrument:            SJI
     Bandpass:              1330.0
-    Obs. Start:            2021-10-01T06:09:24.920
-    Obs. End:              2021-10-01T06:11:44.461
-    Instance Start:        2021-10-01T06:09:25.020
-    Instance End:          2021-10-01T06:11:37.580
-    Total Frames in Obs.:  None
-    IRIS Obs. id:          3683602040
-    IRIS Obs. Description: Very large sparse 16-step raster 15x175 16s   Deep x 0.5 Spatial x 2
+    Obs Date:              2021-10-01T06:09:25.020 -- 2021-10-01T06:11:37.580
+    Total Frames in Obs:   None
+    Obs ID:                3683602040
+    Obs Description:       Very large sparse 16-step raster 15x175 16s   Deep x 0.5 Spatial x 2
     Axis Types:            [('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat', 'time', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM', 'custom:CUSTOM'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat'), ('custom:pos.helioprojective.lon', 'custom:pos.helioprojective.lat')]
     Roll:                  0.000464606
     Cube dimensions:       (20, 548, 555)

examples/offset_sji_sg.py

@@ -0,0 +1,130 @@
+"""
+========================================
+Potential Offset between SJI and SG data
+========================================
+
+In this example, we are showcase a potential offset in WCS coordinates between IRIS SJI and IRIS SG data.
+The offset varies between observations and is not fixed.
+
+The cause is unknown at this time and this has not been cross-checked with SSWIDL.
+So it is possible this is just a bug in the Python library.
+"""
+
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pooch
+
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+from astropy.visualization import time_support
+from astropy.wcs.utils import wcs_to_celestial_frame
+
+from irispy.io import read_files
+
+time_support()
+
+###############################################################################
+# We will start by getting some data from the IRIS archive.
+#
+# In this case, we will use ``pooch`` so to keep this example self-contained
+# but using your browser will also work.
+#
+# Using the url: https://www.lmsal.com/hek/hcr?cmd=view-event&event-id=ivo%3A%2F%2Fsot.lmsal.com%2FVOEvent%23VOEvent_IRIS_20130902_182935_4000005156_2013-09-02T18%3A29%3A352013-09-02T18%3A29%3A35.xml
+# we are after the 2796 Slit-Jaw and the raster sequence.
+
+raster_filename = pooch.retrieve(
+    "https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2013/09/02/20130902_182935_4000005156/iris_l2_20130902_182935_4000005156_raster.tar.gz",
+    known_hash="49444e59cccb5c9b8a41de6fb838de3e43928594685eb0b5f773cbd36a034d86",
+)
+sji_filename = pooch.retrieve(
+    "https://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2013/09/02/20130902_182935_4000005156/iris_l2_20130902_182935_4000005156_SJI_2796_t000_deconvolved.fits.gz",
+    known_hash="19588419fc270aacd296e219d183611ed845a55467e69ad431e68f28ca7e6a53",
+)
+
+###############################################################################
+# Now to open the files using ``irispy``.
+
+raster = read_files(raster_filename)
+sji_2796 = read_files(sji_filename)
+
+###############################################################################
+# Now we will find the closest SJI time to the 56th raster step.
+
+c_ii = raster["C II 1336"][0]
+
+times_SG = c_ii.axis_world_coords("time", wcs=c_ii.extra_coords)
+(time_2796,) = sji_2796.axis_world_coords("time")
+# We picked randomly.
+raster_idx = 55
+time_target = times_SG[0][raster_idx]
+
+time_idx_2796 = np.abs(time_2796 - time_target).argmin()
+time_stamp_2796 = time_2796[time_idx_2796].isot
+print(time_stamp_2796, "\n", time_target, "\n", time_idx_2796)
+
+###############################################################################
+# We will require the axillary data from both files later on/
+
+# The raster file is extracted to the cache pooch via pooch, so we need to do some magic to get it:
+raster_aux_data = fits.getdata(
+    next(iter(Path("~/.cache/pooch/").expanduser().glob("*iris_l2_20130902_182935_4000005156_raster/*fits"))), ext=-2
+)
+sji_aux_data = fits.getdata(sji_filename, ext=-2)
+
+# If you check the headers, you will find an index for POFFXSJI and POFFYSJI:
+# SJI X (spectral direction) shift in pixels
+# SJI Y (spatial direction) shift in pixels
+sji_offset = sji_aux_data[time_idx_2796, 29] * 0.167 * u.arcsec
+sg_offset = raster_aux_data[raster_idx, 45] * 0.167 * u.arcsec
+
+# Also in the SJI AUX data is the pixel location of the slit
+# We will need to get the world location later on.
+sji_slit_location_pixel_x = sji_aux_data[time_idx_2796, 4]
+sji_slit_location_pixel_y = sji_aux_data[time_idx_2796, 5]
+
+###############################################################################
+# We can now get the slit location from the SG cube and apply the offsets to the location.
+
+sji_2796_closest = sji_2796[time_idx_2796]
+sji_2796_frame = wcs_to_celestial_frame(sji_2796_closest.basic_wcs)
+
+# Seems to be a plus one missing?! This seems like a bug. TODO: WORK THIS OUT
+raster_lon_coords = c_ii.axis_world_coords_values("custom:pos.helioprojective.lon")[0][raster_idx + 1]
+raster_lat_coords = c_ii.axis_world_coords_values("custom:pos.helioprojective.lat")[0][raster_idx + 1]
+
+# Now we will get the raster location from its WCS and overlay that on the SJI below.
+slit = SkyCoord(Tx=raster_lon_coords, Ty=raster_lat_coords, frame=sji_2796_frame)
+# In this case, the offsets are negative, so we add them instead.
+slit_with_sg_offset = SkyCoord(Tx=raster_lon_coords + (-1 * sg_offset), Ty=raster_lat_coords, frame=sji_2796_frame)
+slit_with_sji_offset = SkyCoord(Tx=raster_lon_coords + (-1 * sji_offset), Ty=raster_lat_coords, frame=sji_2796_frame)
+
+###############################################################################
+# Now we can visualize the difference.
+
+fig = plt.figure(figsize=(7, 7))
+ax = fig.add_subplot(projection=sji_2796_closest)
+sji_2796_closest.plot(ax)
+
+# This is the pixel location of the slit in the SJI data based on the axillary data
+slit_location_from_sji_aux = sji_2796[time_idx_2796].wcs.pixel_to_world(
+    sji_slit_location_pixel_x, sji_slit_location_pixel_y
+)
+ax.plot_coord(slit_location_from_sji_aux, ".", color="white", label="Slit Pixel Location")
+
+# Now these are the slit locations (some modified by offset values in the axillary data)
+ax.plot_coord(slit, color="white", linestyle="--", linewidth=0.5, label="Slit WCS")
+ax.plot_coord(slit_with_sg_offset, color="red", linestyle="--", linewidth=0.5, label="SG Offset")
+ax.plot_coord(slit_with_sji_offset, color="green", linestyle="--", linewidth=0.5, label="SJI Offset")
+
+ax.legend()
+
+# "Crop" the image without touching the actual data.
+bbox = SkyCoord([140, 180] * u.arcsec, [50, 90] * u.arcsec, frame=sji_2796_frame)
+x_limit, y_limit = sji_2796_closest.wcs.world_to_pixel(bbox)
+ax.set_xlim(int(x_limit[0]), int(x_limit[1]))
+ax.set_ylim(int(y_limit[0]), int(y_limit[1]))
+
+plt.show()

examples/umbral_flashes.py

@@ -26,7 +26,7 @@
 # In this case, we will use ``pooch`` so to keep this example self-contained
 # but using your browser will also work.
 #
-# Using the url: http://www.lmsal.com/solarsoft/irisa/data/level2_compressed/2013/09/02/20130902_163935_4000255147/
+# Using the url: https://www.lmsal.com/hek/hcr?cmd=view-event&event-id=ivo%3A%2F%2Fsot.lmsal.com%2FVOEvent%23VOEvent_IRIS_20130902_163935_4000255147_2013-09-02T16%3A39%3A352013-09-02T16%3A39%3A35.xml
 # we are after the 1400 Slit-Jaw and the raster sequence (~900 MB).
 
 raster_filename = pooch.retrieve(

irispy/io/sji.py

@@ -48,14 +48,12 @@ def _create_gwcs(hdulist: fits.HDUList) -> gwcs.WCS:
         crval_table=crval_table * u.arcsec,
         **kwargs,
     )
-    base_time = Time(hdulist[0].header["STARTOBS"], format="isot", scale="utc")
-    times = hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
-    # We need to account for a non-zero time delta.
-    base_time += times[0]
-    times -= times[0]
+    start_time = Time(hdulist[0].header["DATE_OBS"], format="isot", scale="utc")
+    cadence = hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
+    cadence -= cadence[0]
     temporal = m.Tabular1D(
         np.arange(hdulist[1].data.shape[0]) * u.pix,
-        lookup_table=times,
+        lookup_table=cadence,
         fill_value=np.nan,
         bounds_error=False,
         method="linear",
@@ -64,14 +62,14 @@ def _create_gwcs(hdulist: fits.HDUList) -> gwcs.WCS:
     celestial_frame = cf.CelestialFrame(
         axes_order=(0, 1),
         unit=(u.arcsec, u.arcsec),
-        reference_frame=Helioprojective(observer="earth", obstime=base_time),
+        reference_frame=Helioprojective(observer="earth", obstime=start_time),
         axis_physical_types=[
             "custom:pos.helioprojective.lon",
             "custom:pos.helioprojective.lat",
         ],
         axes_names=("Longitude", "Latitude"),
     )
-    temporal_frame = cf.TemporalFrame(Time(base_time), unit=(u.s,), axes_order=(2,), axes_names=("Time (UTC)",))
+    temporal_frame = cf.TemporalFrame(start_time, unit=(u.s,), axes_order=(2,), axes_names=("Time (UTC)",))
     output_frame = cf.CompositeFrame([celestial_frame, temporal_frame])
     input_frame = cf.CoordinateFrame(
         axes_order=(0, 1, 2),
@@ -90,11 +88,10 @@ def _create_wcs(hdulist):
     This has been set to have an Earth Observer at the time of the observation.
     """
     wcses = []
-    base_time = Time(hdulist[0].header["STARTOBS"], format="isot", scale="utc")
-    times = hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
-    # We need to account for a non-zero time delta.
-    base_time += times[0]
-    times -= times[0]
+    obs_times = (
+        Time(hdulist[0].header["STARTOBS"], format="isot", scale="utc")
+        + hdulist[1].data[:, hdulist[1].header["TIME"]] * u.s
+    )
     xcenix_idx = hdulist[1].header["XCENIX"]
     ycenix_idx = hdulist[1].header["YCENIX"]
     pc1_1ix = hdulist[1].header["PC1_1IX"]
@@ -122,12 +119,12 @@ def _create_wcs(hdulist):
             nonzero_vals = array[nonzero_idx]
             array[zero_idx] = np.interp(zero_idx, nonzero_idx, nonzero_vals)
     for i in range(hdulist[0].header["NAXIS3"]):
-        location = get_body_heliographic_stonyhurst("Earth", (base_time + times[i]).isot)
+        location = get_body_heliographic_stonyhurst("Earth", (obs_times[i]).isot)
         observer = Helioprojective(
             xcenix_values[i] * u.arcsec,
             ycenix_values[i] * u.arcsec,
             observer=location,
-            obstime=base_time + times[i],
+            obstime=obs_times[i],
         )
         rotation_matrix = np.asanyarray(
             [

irispy/io/spectrograph.py

@@ -92,7 +92,7 @@ def read_spectrograph_lvl2(
             window_fits_indices = np.nonzero(np.isin(windows_in_obs, spectral_windows))[0] + 1
         data_dict = {window_name: [] for window_name in spectral_windows_req}
         # No observer information in the header, so we just assume its at Earth.
-        base_time = Time(hdulist[0].header["STARTOBS"])
+        base_time = Time(hdulist[0].header["DATE_OBS"])
         location = get_body_heliographic_stonyhurst("Earth", (base_time).isot)
         observer = Helioprojective(
             hdulist[0].header["XCEN"] * u.arcsec,

irispy/sji.py

@@ -95,13 +95,10 @@ def __str__(self) -> str:
             Observatory:           {self.meta.get("TELESCOP", "IRIS")}
             Instrument:            {self.meta.get("INSTRUME")}
             Bandpass:              {self.meta.get("TWAVE1")}
-            Obs. Start:            {self.meta.get("STARTOBS")}
-            Obs. End:              {self.meta.get("ENDOBS")}
-            Instance Start:        {instance_start}
-            Instance End:          {instance_end}
-            Total Frames in Obs.:  {self.meta.get("NBFRAMES")}
-            IRIS Obs. id:          {self.meta.get("OBSID")}
-            IRIS Obs. Description: {self.meta.get("OBS_DESC")}
+            Obs Date:              {instance_start} -- {instance_end}
+            Total Frames in Obs:   {self.meta.get("NBFRAMES")}
+            Obs ID:                {self.meta.get("OBSID")}
+            Obs Description:       {self.meta.get("OBS_DESC")}
             Axis Types:            {self.array_axis_physical_types}
             Roll:                  {self.meta.get("SAT_ROT")}
             Cube dimensions:       {self.shape}

irispy/spectrograph.py

@@ -81,10 +81,9 @@ def __str__(self) -> str:
             f"""
             SpectrogramCube
             ---------------
-            OBS ID:             {self.meta.get("OBSID")}
-            OBS Description:    {self.meta.get("OBS_DESC")}
-            OBS period:         {self.meta.get("STARTOBS")} -- {self.meta.get("ENDOBS")}
-            Spectrogram period: {instance_start} -- {instance_end}
+            Obs ID:             {self.meta.get("OBSID")}
+            Obs Description:    {self.meta.get("OBS_DESC")}
+            Obs Date:           {instance_start} -- {instance_end}
             Data shape:         {self.shape}
             Axis Types:         {self.array_axis_physical_types}
             Roll:               {self.meta.get("SAT_ROT")}

irispy/utils/constants.py

@@ -36,12 +36,20 @@
 RADIANCE_UNIT = u.erg / u.cm**2 / u.s / u.steradian / u.Angstrom
 SLIT_WIDTH = 0.33 * u.arcsec
 SPECTRAL_BAND = {
+    "1336": "FUV",
     "1343": "FUV",
+    "1349": "FUV",
+    "1352": "FUV",
+    "1356": "FUV",
+    "1394": "FUV",
     "1400": "FUV",
+    "1403": "FUV",
     "2786": "NUV",
+    "2796": "NUV",
     "2814": "NUV",
     "2826": "NUV",
     "2830": "NUV",
+    "2831": "NUV",
     "2832": "NUV",
     "C II 1336": "FUV",
     "Cl I 1352": "FUV",

irispy/utils/wobble.py

@@ -68,7 +68,7 @@ def generate_wobble_movie(
 
     2832 is considered the best wavelength to use for wobble movies.
 
-    Timestamps take the main header cadence and add that to the "DATEOBS".
+    Timestamps take the main header cadence and add that to the "DATE_OBS".
     They do not use the information in the AUX array.
     """
     # Avoid circular imports
@@ -90,10 +90,10 @@ def generate_wobble_movie(
         cadence_sample = max(1, np.floor(wobble_cadence / cadence))
         if timestamp:
             timestamps = [
-                parse_time(header["STARTOBS"]) + TimeDelta(cadence, format="sec") * i for i in range(numframes)
+                parse_time(header["DATE_OBS"]) + TimeDelta(cadence, format="sec") * i for i in range(numframes)
             ]
         else:
-            timestamps = [parse_time(header["STARTOBS"])]
+            timestamps = [parse_time(header["DATE_OBS"])]
         # Trim down to only that part of the movie that contains data in all frames
         if trim:
             # TODO: improve this, it trims a bit but not fully

ruff.toml

@@ -55,6 +55,7 @@ lint.extend-ignore = [
     "ERA001", # Commented out code
     "INP001", # Implicit namespace package
     "T201", # Use print
+    "N816", # Variable `times_SG` in global scope should not be mixedCase"
 ]
 "docs/conf.py" = [
     "INP001", # conf.py is part of an implicit namespace package

tox.ini

@@ -1,7 +1,7 @@
 [tox]
 min_version = 4.0
 envlist =
-    py{311,312,313}{,-online,-figure}
+    py{311,312,313,314}{,-online,-figure}
     py313-devdeps
     py311-oldestdeps
     codestyle