adding xrt_teem

.gitignore

@@ -141,3 +141,6 @@ version.py
 #ignore mac files
 .DS_store
 .DS_Store
+
+# local directory
+idl_routines/

docs/api_static/xrtpy.response.xrt_teem.rst

@@ -0,0 +1,8 @@
+:orphan:
+
+`xrtpy.response.xrt_teem`
+=========================
+
+.. currentmodule:: xrtpy.response.xrt_teem
+
+.. automodapi:: xrtpy.response.xrt_teem

docs/examples.rst

@@ -23,3 +23,11 @@ Computing response functions
   :glob:
 
   notebooks/computing_functions/*
+
+Data analysis
+-------------
+
+.. nbgallery::
+   :glob:
+
+   Using_xrt_teem

docs/getting_started.rst

@@ -34,6 +34,10 @@ The default model assumes coronal abundances :cite:t:`feldman:1992`.
 .. note::
    XRTpy has future plans to accept other plasma emission spectra models.
 
+Deriving Temperature and Emission Measure for a Pair of Images
+--------------------------------------------------------------
+XRTpy provides a routine, xrt_teem, that employs the objects listed above to derive the temperature and emission measure in for a given pair of images using the filter ratio method. This uses the same methods as in the SolarSoft IDL routine of the same name.
+
 
 X-Ray Filter Channel
 *********************

docs/notebooks/data_analysis/Using_xrt_teem.ipynb

@@ -0,0 +1,191 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "642ac1b3",
+   "metadata": {},
+   "source": [
+    "# Using xrt_teem to analyze XRT data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "22b53fce",
+   "metadata": {},
+   "source": [
+    "This example demonstrates how to use xrt_teem to calculate the temperature and emission measure in an image using the filter ratio method."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f5f04496",
+   "metadata": {},
+   "source": [
+    "First we need to import xrt_teem."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "58758bf1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xrtpy.response.xrt_teem import xrt_teem"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8aee2bae",
+   "metadata": {},
+   "source": [
+    "As an example we will use the test data included in XRTpy, though data with the right characteristics in the XRT archive could also be used. It's necessary to use two images that are the same size and use different filters. To get good results the images should have been taken close in time as well, ideally adjacent in time. Note that not all filter ratios produce good results. \n",
+    "\n",
+    "This data was generated using the IDL routine xrt_prep.pro from SolarSoft and is unnormalized. Data in the Level 1 archive are normalized, which is also okay to use, though the IDL routine xrt_teem.pro did not allow that. For normalized data the image data is multiplied by the exposure time before analysis. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da83268f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pkg_resources\n",
+    "import sunpy.map\n",
+    "\n",
+    "from pathlib import Path\n",
+    "\n",
+    "directory = pkg_resources.resource_filename(\n",
+    "    \"xrtpy\", \"response/tests/data/xrt_teem_testing_files\"\n",
+    ")\n",
+    "data_files = sorted(Path(directory).glob(\"L1_XRT20110128_*.*.fits\"))\n",
+    "file1 = data_files[1]\n",
+    "file2 = data_files[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "35dd0c83",
+   "metadata": {},
+   "source": [
+    "xrt_teem uses SunPy maps as input. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "48ac4220",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "map1 = sunpy.map.Map(file1)\n",
+    "map2 = sunpy.map.Map(file2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "905f4e11",
+   "metadata": {},
+   "source": [
+    "xrt_teem has several options, mirroring the IDL routine in SolarSoft in most respects. A simple call with no extra parameters calculates the temperature and (volume) emission measure for the two images without any binning or masking of the data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ddc9862b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "T_e, EM, Terr, EMerr = xrt_teem(map1, map2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d0ef85f3",
+   "metadata": {},
+   "source": [
+    "The outputs of xrt_teem are also SunPy maps. T_e is the electron temperature, EM is the volume emission measure, Terr is a measure of the uncertainties in the temperature determined for each pixel and EMerr is the same for the emission measure. Each map has data and associated metadata. To examine the results one can use matplotlib and sunpy:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3ae1f53b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from sunpy.coordinates.sun import angular_radius, B0\n",
+    "from sunpy.map import Map\n",
+    "\n",
+    "# To avoid error messages from sunpy we add metadata to the header:\n",
+    "rsun_ref = 6.95700e08\n",
+    "hdr1 = map1.meta\n",
+    "rsun_obs = angular_radius(hdr1[\"DATE_OBS\"]).value\n",
+    "dsun = rsun_ref / np.sin(rsun_obs * np.pi / 6.48e5)\n",
+    "solarb0 = B0(hdr1[\"DATE_OBS\"]).value\n",
+    "hdr1[\"DSUN_OBS\"] = dsun\n",
+    "hdr1[\"RSUN_REF\"] = rsun_ref\n",
+    "hdr1[\"RSUN_OBS\"] = rsun_obs\n",
+    "hdr1[\"SOLAR_B0\"] = solarb0\n",
+    "\n",
+    "fig = plt.figure()\n",
+    "# We could create a plot simply by doing T_e.plot(), but here we choose to make a linear plot of T_e\n",
+    "m = Map((10.0**T_e.data, T_e.meta))\n",
+    "m.plot(title=\"Derived Temperature\", vmin=2.0e6, vmax=1.2e7, cmap=\"turbo\")\n",
+    "m.draw_limb()\n",
+    "m.draw_grid(linewidth=2)\n",
+    "cb = plt.colorbar(label=\"T (K)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4e069966",
+   "metadata": {},
+   "source": [
+    "See the xrt_teem.py code for more information. Among the options are verbose output, binning the data by an integer factor (to increase the signal to noise), specifying a temperature range to examine, providing a mask for excluding regions of the images from the analysis, and setting error thresholds on the temperature and photon noise that differ from the default values."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "94df8fee",
+   "metadata": {},
+   "source": [
+    "These data were analyzed by Guidoni et al. (2015, ApJ 800, 54). See also Narukage et al. (2014, Solar Phys. 289, 1029). The filters for these two images are Ti poly and Be thin. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2807a18e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

setup.cfg

@@ -120,6 +120,7 @@ extend-ignore =
     W605,
     RST210,
     RST213,
+    RST305,
     RST306
 
 exclude =

xrtpy/response/tests/test_xrt_teem.py

@@ -0,0 +1,141 @@
+import numpy as np
+import pkg_resources
+import pytest
+import sunpy.map
+
+from astropy.io import fits
+from pathlib import Path
+from scipy.io import readsav
+
+from xrtpy.response.xrt_teem import xrt_teem
+
+
+def get_observed_data():
+    directory = pkg_resources.resource_filename(
+        "xrtpy", "response/tests/data/xrt_teem_testing_files"
+    )
+    data_files = sorted(Path(directory).glob("L1_XRT20110128_*.*.fits"))
+
+    return data_files
+
+
+def get_IDL_results_data():
+    directory = pkg_resources.resource_filename(
+        "xrtpy", "response/tests/data/xrt_teem_testing_files"
+    )
+    results_files = sorted(Path(directory).glob("IDL_results_*.sav"))
+
+    return results_files
+
+
+def rebin_image(data, binfac=1):
+    """
+    Given a data array and a binning factor return the data array rebinned by
+    the binning factor.
+    """
+
+    s = data.shape
+    ns = (s[0] // binfac, s[1] // binfac)
+    rbs = (ns[0], binfac, ns[1], binfac)
+    # sums the data in binfac x binfac sized regions
+    drbin = data.reshape(rbs).mean(-1).mean(1)
+    # for a boolean mask, this makes a pixel masked if any of the summed
+    # pixels is masked. If we want to mask only if all the pixels are masked
+    # then we could use prod in place of sum above
+
+    if data.dtype == bool:
+        # need to convert back to bool after summing
+        drbin = drbin.astype(bool)
+    return drbin
+
+
+def test_standard_case():
+    """
+    Test case with all default values:
+    no binning, no masking, no temperature range, standard thresholds
+    """
+
+    data_files = get_observed_data()
+    # it turns out that the IDL test data was generated with the inverse order
+    # of the files
+    file1 = data_files[1]
+    file2 = data_files[0]
+    map1 = sunpy.map.Map(file1)
+    map2 = sunpy.map.Map(file2)
+
+    T_e, EM, Terr, EMerr = xrt_teem(map1, map2)
+
+    testdata = get_IDL_results_data()
+
+    # This is needed because there are multiple test data sets, though so far
+    # only a test written for the standard case
+    fnames = [td.name for td in testdata]
+    idata1 = fnames.index("IDL_results_bin1.sav")
+    testdata1 = testdata[idata1]
+
+    idldata = readsav(testdata1)
+    goodT = (T_e.data > 0.0) & (idldata.te > 0.0)
+    goodE = (EM.data > 0.0) & (idldata.em > 0.0)
+    assert np.allclose(
+        10.0 ** T_e.data[goodT], 10.0 ** idldata.te[goodT], atol=2.0e5, rtol=0.02
+    )
+    assert np.allclose(
+        10.0 ** EM.data[goodE], 10.0 ** idldata.em[goodE], atol=4.0e44, rtol=0.03
+    )
+    assert np.allclose(
+        10.0 ** Terr.data[goodT], 10.0 ** idldata.et[goodT], atol=1.0e4, rtol=0.08
+    )
+    assert np.allclose(
+        10.0 ** EMerr.data[goodE], 10.0 ** idldata.ee[goodE], atol=4.0e43, rtol=0.02
+    )
+
+
+def test_binning_case():
+    """
+    Test case with following parameters:
+    binning by a factor of 2
+    no masking
+    no temperature range
+    standard thresholds
+    """
+
+    data_files = get_observed_data()
+    # it turns out that the IDL test data was generated with the inverse order
+    # of the files
+    file1 = data_files[1]
+    file2 = data_files[0]
+    map1 = sunpy.map.Map(file1)
+    map2 = sunpy.map.Map(file2)
+
+    T_e, EM, Terr, EMerr = xrt_teem(map1, map2, binfac=2)
+
+    testdata = get_IDL_results_data()
+
+    # This is needed because there are multiple test data sets
+    fnames = [td.name for td in testdata]
+    idata1 = fnames.index("IDL_results_bin2.sav")
+    testdata1 = testdata[idata1]
+
+    idldata = readsav(testdata1)
+    idlTe = rebin_image(idldata.te, 2)
+    idlEM = rebin_image(idldata.em, 2)
+    idlTerr = rebin_image(idldata.et, 2)
+    idlEMerr = rebin_image(idldata.ee, 2)
+    goodT = (T_e.data > 0.0) & (idlTe > 0.0)
+    goodE = (EM.data > 0.0) & (idlEM > 0.0)
+
+    delta = 10.0 ** T_e.data[goodT] - 10.0 ** idlTe[goodT]
+    x = 10.0 ** idlTe[goodT]
+
+    assert np.allclose(
+        10.0 ** T_e.data[goodT], 10.0 ** idlTe[goodT], atol=2.0e5, rtol=0.02
+    )
+    assert np.allclose(
+        10.0 ** EM.data[goodE], 10.0 ** idlEM[goodE], atol=1.0e44, rtol=0.05
+    )
+    assert np.allclose(
+        10.0 ** Terr.data[goodT], 10.0 ** idlTerr[goodT], atol=1.0e4, rtol=0.1
+    )
+    assert np.allclose(
+        10.0 ** EMerr.data[goodE], 10.0 ** idlEMerr[goodE], atol=2.0e43, rtol=0.03
+    )