From 8ab75f687dfbcfd7a3f3225245aade8b566ed0e1 Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Fri, 12 Dec 2025 11:26:08 -0800 Subject: [PATCH] Added Offset example and tweaked times --- docs/guide.rst | 11 ++-- examples/offset_sji_sg.py | 130 +++++++++++++++++++++++++++++++++++++ examples/umbral_flashes.py | 2 +- irispy/io/sji.py | 27 ++++---- irispy/io/spectrograph.py | 2 +- irispy/sji.py | 11 ++-- irispy/spectrograph.py | 7 +- irispy/utils/constants.py | 8 +++ irispy/utils/wobble.py | 6 +- ruff.toml | 1 + tox.ini | 2 +- 11 files changed, 168 insertions(+), 39 deletions(-) create mode 100644 examples/offset_sji_sg.py diff --git a/docs/guide.rst b/docs/guide.rst index 13d2d60..b664e6c 100644 --- a/docs/guide.rst +++ b/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) diff --git a/examples/offset_sji_sg.py b/examples/offset_sji_sg.py new file mode 100644 index 0000000..ed3d969 --- /dev/null +++ b/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() diff --git a/examples/umbral_flashes.py b/examples/umbral_flashes.py index 0d8df21..f456a0e 100644 --- a/examples/umbral_flashes.py +++ b/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( diff --git a/irispy/io/sji.py b/irispy/io/sji.py index f82a487..37b7dbc 100644 --- a/irispy/io/sji.py +++ b/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( [ diff --git a/irispy/io/spectrograph.py b/irispy/io/spectrograph.py index 010af40..ccb33c2 100644 --- a/irispy/io/spectrograph.py +++ b/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, diff --git a/irispy/sji.py b/irispy/sji.py index 185c32c..ddb4f88 100644 --- a/irispy/sji.py +++ b/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} diff --git a/irispy/spectrograph.py b/irispy/spectrograph.py index 34cc7aa..29c8fb0 100644 --- a/irispy/spectrograph.py +++ b/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")} diff --git a/irispy/utils/constants.py b/irispy/utils/constants.py index 0d3159c..fc9970d 100644 --- a/irispy/utils/constants.py +++ b/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", diff --git a/irispy/utils/wobble.py b/irispy/utils/wobble.py index d1139ba..384648a 100644 --- a/irispy/utils/wobble.py +++ b/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 diff --git a/ruff.toml b/ruff.toml index 3fc45a7..41db6de 100644 --- a/ruff.toml +++ b/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 diff --git a/tox.ini b/tox.ini index d32f785..cded3ef 100644 --- a/tox.ini +++ b/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