From 391b935bd8ca51517cb677ceb641a55bf5189460 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 5 Nov 2024 19:05:07 -0700 Subject: [PATCH] Added `iris.sg.SpectrographObservation.radiance` property which converts from DN to radiometric units. (#21) --- iris/sg/__init__.py | 8 ++++++- iris/sg/_effective_area.py | 14 ++++++++++++ iris/sg/_spectrograph.py | 40 +++++++++++++++++++++++++++++++++++ iris/sg/_spectrograph_test.py | 8 +++++++ 4 files changed, 69 insertions(+), 1 deletion(-) diff --git a/iris/sg/__init__.py b/iris/sg/__init__.py index 3eb0ad3..737d460 100644 --- a/iris/sg/__init__.py +++ b/iris/sg/__init__.py @@ -3,11 +3,17 @@ """ from . import background -from ._effective_area import effective_area +from ._effective_area import ( + gain, + width_slit, + effective_area, +) from ._spectrograph import SpectrographObservation __all__ = [ "background", + "gain", + "width_slit", "effective_area", "SpectrographObservation", ] diff --git a/iris/sg/_effective_area.py b/iris/sg/_effective_area.py index 057e680..388f40a 100644 --- a/iris/sg/_effective_area.py +++ b/iris/sg/_effective_area.py @@ -4,9 +4,23 @@ import iris __all__ = [ + "gain", + "width_slit", "effective_area", ] +gain = 12 * u.photon / u.DN +""" +The conversion factor between counts and photons measured by the sensor +:cite:p:`Wulser2018`. +""" + + +width_slit = 1 / 3 * u.arcsec +""" +The angular subtent of the spectrographic slit. +""" + def effective_area(wavelength: u.Quantity | na.AbstractScalar) -> na.AbstractScalar: """ diff --git a/iris/sg/_spectrograph.py b/iris/sg/_spectrograph.py index 3df043f..628da79 100644 --- a/iris/sg/_spectrograph.py +++ b/iris/sg/_spectrograph.py @@ -1,7 +1,9 @@ +from typing_extensions import Self import dataclasses import pathlib import numpy as np import astropy.units as u +import astropy.constants import astropy.time import astropy.wcs import astropy.io.fits @@ -432,3 +434,41 @@ def empty( axis_detector_x=axis_detector_x, axis_detector_y=axis_detector_y, ) + + @property + def radiance(self) -> Self: + """ + Convert to radiometric units using :func:`iris.sg.effective_area`. + """ + + wavelength = self.inputs.wavelength + + lower = {self.axis_wavelength: slice(None, ~0)} + upper = {self.axis_wavelength: slice(+1, None)} + wavelength = (wavelength[lower] + wavelength[upper]) / 2 + + energy = astropy.constants.h * astropy.constants.c / wavelength / u.ph + + gain = iris.sg.gain + + area_eff = iris.sg.effective_area(wavelength) + + pix_xy = np.diff(self.inputs.position, axis=self.axis_detector_y).length + pix_xy = pix_xy.mean(self.axis_detector_x) + + pix_lambda = np.diff(self.inputs.wavelength, axis=self.axis_wavelength) + + t_exp = self.timedelta + + w_slit = iris.sg.width_slit + + factor = energy * gain / (area_eff * pix_xy * pix_lambda * t_exp * w_slit) + + outputs = self.outputs * factor + + outputs = outputs.to(u.erg / (u.cm**2 * u.sr * u.nm * u.s)) + + return dataclasses.replace( + self, + outputs=outputs, + ) diff --git a/iris/sg/_spectrograph_test.py b/iris/sg/_spectrograph_test.py index 56fc5dd..79b1fc3 100644 --- a/iris/sg/_spectrograph_test.py +++ b/iris/sg/_spectrograph_test.py @@ -1,4 +1,6 @@ import pytest +import numpy as np +import astropy.units as u import astropy.time import iris @@ -25,3 +27,9 @@ def test_axis_detector_x(self, array: iris.sg.SpectrographObservation): def test_axis_detector_y(self, array: iris.sg.SpectrographObservation): assert isinstance(array.axis_detector_y, str) + + def test_radiance(self, array: iris.sg.SpectrographObservation): + result = array.radiance + assert isinstance(result, iris.sg.SpectrographObservation) + assert np.all(result.inputs == array.inputs) + assert np.nansum(result.outputs) > 0 * u.erg / (u.cm**2 * u.sr * u.s * u.nm)