changeset 0:c80db1ec6611 draft default tip

planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit 5b1251e9c9cdff7f825bfcfebdfdb12714c13d8f
author astroteam
date Wed, 13 Aug 2025 19:10:53 +0000
parents
children
files detectgrb.py grb_detection_astro_tool.xml
diffstat 2 files changed, 307 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/detectgrb.py	Wed Aug 13 19:10:53 2025 +0000
@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+#!/usr/bin/env python
+
+# This script is generated with nb2galaxy
+
+# flake8: noqa
+
+import json
+import os
+import shutil
+
+from oda_api.json import CustomJSONEncoder
+
+T1 = "2023-01-16T04:53:33.9"  # http://odahub.io/ontology#StartTime
+T2 = "2023-01-16T04:55:33.9"  # http://odahub.io/ontology#EndTime
+detection_time_scales = "1,10"
+lc_time_scale = 0.1  # https://odahub.io/ontology#TimeIntervalSeconds
+background_age = 10  # oda:TimeIntervalSeconds
+min_sn = 5  # https://odahub.io/ontology#SignalToNoiseRatio
+
+_galaxy_wd = os.getcwd()
+
+with open("inputs.json", "r") as fd:
+    inp_dic = json.load(fd)
+if "C_data_product_" in inp_dic.keys():
+    inp_pdic = inp_dic["C_data_product_"]
+else:
+    inp_pdic = inp_dic
+T1 = str(inp_pdic["T1"])
+T2 = str(inp_pdic["T2"])
+detection_time_scales = str(inp_pdic["detection_time_scales"])
+lc_time_scale = float(inp_pdic["lc_time_scale"])
+background_age = float(inp_pdic["background_age"])
+min_sn = int(inp_pdic["min_sn"])
+
+import numpy as np
+from astropy.time import Time
+from matplotlib import pylab as plt
+from oda_api.api import DispatcherAPI
+from oda_api.data_products import LightCurveDataProduct, PictureProduct
+
+disp = DispatcherAPI(
+    url="https://www.astro.unige.ch/mmoda//dispatch-data", instrument="mock"
+)
+
+par_dict = {
+    "T1": T1,
+    "T2": T2,
+    "T_format": "isot",
+    "instrument": "spi_acs",
+    "product": "spi_acs_lc",
+    "product_type": "Real",
+    "time_bin": lc_time_scale,
+    "time_bin_format": "sec",
+}
+
+data_collection = disp.get_product(**par_dict)
+
+lc = data_collection.spi_acs_lc_0_query.data_unit[1].data
+
+lc["TIME"] = (
+    lc["TIME"] / 24 / 3600 + (Time(T1).mjd + Time(T2).mjd) / 2.0
+)  # TODO: more accurately
+
+def rebin(x, n):
+    N = int(len(x) / n)
+    return np.array(x[: N * n]).reshape((N, n)).sum(1)
+
+lc = LightCurveDataProduct.from_arrays(
+    times=Time(lc["TIME"], format="mjd"), fluxes=lc["RATE"], errors=lc["ERROR"]
+)
+
+obj_results = [lc.encode()]
+
+T0_ijd = lc.data_unit[1].data["TIME"][np.argmax(lc.data_unit[1].data["FLUX"])]
+
+ijd2plot = lambda x: (x - T0_ijd) * 24 * 3600
+
+m_bkg = ijd2plot(lc.data_unit[1].data["TIME"]) < background_age
+np.sum(m_bkg)
+bkg = np.mean(lc.data_unit[1].data["FLUX"][m_bkg])
+bkg
+
+from matplotlib import pylab as plt
+
+plt.figure(figsize=(10, 6))
+
+plt.errorbar(
+    ijd2plot(lc.data_unit[1].data["TIME"]),
+    lc.data_unit[1].data["FLUX"],
+    lc.data_unit[1].data["ERROR"],
+    label=f"{lc_time_scale} s",
+    alpha=0.5,
+)
+
+best_detection = {"sn": None, "timescale": None, "t_max_sn": None}
+
+for detection_time_scale in detection_time_scales.split(","):
+    n = int(float(detection_time_scale.strip()) / lc_time_scale)
+
+    r_t = rebin(lc.data_unit[1].data["TIME"], n) / n
+    r_c = rebin(lc.data_unit[1].data["FLUX"], n) / n
+    r_ce = rebin(lc.data_unit[1].data["ERROR"] ** 2, n) ** 0.5 / n
+
+    sn = (r_c - bkg) / r_ce
+    imax = sn.argmax()
+
+    if best_detection["sn"] is None or best_detection["sn"] < sn[imax]:
+        best_detection["sn"] = sn[imax]
+        best_detection["timescale"] = detection_time_scale
+        best_detection["t_max_sn"] = Time(r_t[imax], format="mjd").isot
+
+    plt.errorbar(
+        ijd2plot(r_t),
+        r_c,
+        r_ce,
+        label=f"{detection_time_scale} s, S/N = {sn[imax]:.1f}",
+    )
+
+plt.axhline(bkg, lw=3, ls=":", c="k")
+
+plt.grid()
+plt.legend()
+
+plt.ylabel("Flux")
+plt.xlabel("Time, MJD")
+
+plt.title(best_detection)
+
+plt.savefig("annotated-lc.png")
+
+Time(r_t, format="mjd").isot
+
+bin_image = PictureProduct.from_file("annotated-lc.png")
+
+detection_note = str(dict(best_detection=best_detection))
+
+import numpy as np
+from astropy import wcs
+from astropy.io import fits
+from matplotlib import pylab as plt
+
+x, y = np.meshgrid(np.linspace(0, 10, 500), np.linspace(0, 10, 500))
+
+image = np.exp(-((x - 5) ** 2 + (y - 5) ** 2))
+
+pixsize = 0.1
+Npix = 500
+cdec = 1
+
+RA = 0
+DEC = 0
+
+# Create a new WCS object.  The number of axes must be set
+# from the start
+w = wcs.WCS(naxis=2)
+
+# Set up an "Airy's zenithal" projection
+# Vector properties may be set with Python lists, or Numpy arrays
+w.wcs.crpix = [image.shape[0] / 2.0, image.shape[1] / 2.0]
+w.wcs.cdelt = np.array([pixsize / cdec, pixsize])
+w.wcs.crval = [RA, DEC]
+w.wcs.ctype = ["RA---CAR", "DEC--CAR"]
+w.wcs.set_pv([(1, 1, 45.0)])
+
+# Now, write out the WCS object as a FITS header
+header = w.to_header()
+
+# header is an astropy.io.fits.Header object.  We can use it to create a new
+# PrimaryHDU and write it to a file.
+hdu = fits.PrimaryHDU(image, header=header)
+hdu.writeto("Image.fits", overwrite=True)
+hdu = fits.open("Image.fits")
+im = hdu[0].data
+from astropy.wcs import WCS
+
+wcs = WCS(hdu[0].header)
+plt.subplot(projection=wcs)
+plt.imshow(im, origin="lower")
+plt.grid(color="white", ls="solid")
+plt.xlabel("RA")
+plt.ylabel("Dec")
+
+from oda_api.data_products import ImageDataProduct
+
+fits_image = ImageDataProduct.from_fits_file("Image.fits")
+
+lc = lc  # http://odahub.io/ontology#LightCurve
+
+detection_comment = detection_note  # http://odahub.io/ontology#ODATextProduct
+image_output = bin_image  # http://odahub.io/ontology#ODAPictureProduct
+image_fits = fits_image  # http://odahub.io/ontology#Image
+
+# output gathering
+_galaxy_meta_data = {}
+_oda_outs = []
+_oda_outs.append(("out_detectgrb_lc", "lc_galaxy.output", lc))
+_oda_outs.append(
+    (
+        "out_detectgrb_detection_comment",
+        "detection_comment_galaxy.output",
+        detection_comment,
+    )
+)
+_oda_outs.append(
+    ("out_detectgrb_image_output", "image_output_galaxy.output", image_output)
+)
+_oda_outs.append(
+    ("out_detectgrb_image_fits", "image_fits_galaxy.output", image_fits)
+)
+
+for _outn, _outfn, _outv in _oda_outs:
+    _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
+    if isinstance(_outv, str) and os.path.isfile(_outv):
+        shutil.move(_outv, _galaxy_outfile_name)
+        _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
+    elif getattr(_outv, "write_fits_file", None):
+        _outv.write_fits_file(_galaxy_outfile_name)
+        _galaxy_meta_data[_outn] = {"ext": "fits"}
+    elif getattr(_outv, "write_file", None):
+        _outv.write_file(_galaxy_outfile_name)
+        _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
+    else:
+        with open(_galaxy_outfile_name, "w") as fd:
+            json.dump(_outv, fd, cls=CustomJSONEncoder)
+        _galaxy_meta_data[_outn] = {"ext": "json"}
+
+with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
+    json.dump(_galaxy_meta_data, fd)
+print("*** Job finished successfully ***")
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/grb_detection_astro_tool.xml	Wed Aug 13 19:10:53 2025 +0000
@@ -0,0 +1,75 @@
+<tool id="grb_detection_astro_tool" name="GRB detection" version="0.0.1+galaxy0" profile="24.0">
+  <requirements>
+    <requirement type="package" version="1.26.4">numpy</requirement>
+    <requirement type="package" version="6.1.4">astropy</requirement>
+    <requirement type="package" version="3.9.4">matplotlib</requirement>
+    <requirement type="package" version="1.2.41">oda-api</requirement>
+    <requirement type="package" version="3.12.11">python</requirement>
+    <requirement type="package" version="9.4.0">ipython</requirement>
+    <requirement type="package" version="7.16.6">nbconvert</requirement>
+    <requirement type="package" version="5.10.4">nbformat</requirement>
+    <requirement type="package" version="0.3.4">gcn-kafka</requirement>
+    <requirement type="package" version="0.14.2">xmltodict</requirement>
+    <requirement type="package" version="0.11.1">hop-client</requirement>
+    <requirement type="package" version="2.4.0">ligo.skymap</requirement>
+  </requirements>
+  <command detect_errors="exit_code">python '$__tool_directory__/detectgrb.py'</command>
+  <environment_variables>
+    <environment_variable name="BASEDIR">$__tool_directory__</environment_variable>
+    <environment_variable name="GALAXY_TOOL_DIR">$__tool_directory__</environment_variable>
+  </environment_variables>
+  <configfiles>
+    <inputs name="inputs" filename="inputs.json" data_style="paths" />
+  </configfiles>
+  <inputs>
+    <param name="T1" type="text" value="2023-01-16T04:53:33.9" label="T1" optional="false" />
+    <param name="T2" type="text" value="2023-01-16T04:55:33.9" label="T2" optional="false" />
+    <param name="detection_time_scales" type="text" value="1,10" label="detection_time_scales" optional="false" />
+    <param name="lc_time_scale" type="float" value="0.1" label="lc_time_scale" optional="false" />
+    <param name="background_age" type="float" value="10" label="background_age (unit: s)" optional="false" />
+    <param name="min_sn" type="integer" value="5" label="min_sn" optional="false" />
+  </inputs>
+  <outputs>
+    <data label="${tool.name} -&gt; detectgrb lc" name="out_detectgrb_lc" format="auto" from_work_dir="lc_galaxy.output" />
+    <data label="${tool.name} -&gt; detectgrb detection_comment" name="out_detectgrb_detection_comment" format="auto" from_work_dir="detection_comment_galaxy.output" />
+    <data label="${tool.name} -&gt; detectgrb image_output" name="out_detectgrb_image_output" format="auto" from_work_dir="image_output_galaxy.output" />
+    <data label="${tool.name} -&gt; detectgrb image_fits" name="out_detectgrb_image_fits" format="auto" from_work_dir="image_fits_galaxy.output" />
+  </outputs>
+  <tests>
+    <test expect_num_outputs="4">
+      <param name="T1" value="2023-01-16T04:53:33.9" />
+      <param name="T2" value="2023-01-16T04:55:33.9" />
+      <param name="detection_time_scales" value="1,10" />
+      <param name="lc_time_scale" value="0.1" />
+      <param name="background_age" value="10" />
+      <param name="min_sn" value="5" />
+      <assert_stdout>
+        <has_text text="*** Job finished successfully ***" />
+      </assert_stdout>
+    </test>
+  </tests>
+  <help>This service provides simple detection process for a gamma-ray burst in
+SPI-ACS data.
+
+It queries MMODA platform to get SPI-ACS lightcurve in the interval
+between ``T1`` and ``T2`` and detects GRB as a peak above the background
+(calculated in the window of the ``background_age`` width) with given
+minimal S/N ratio. The detection is performed in the rescaled lightcurve
+(scales given as comma-separated list).
+</help>
+  <citations>
+    <citation type="bibtex">@article{Neronov_2021,
+			title = {Online data analysis system of the INTEGRAL telescope},
+			volume = {651},
+			ISSN = {1432-0746},
+			url = {http://dx.doi.org/10.1051/0004-6361/202037850},
+			DOI = {10.1051/0004-6361/202037850},
+			journal = {Astronomy &amp;amp; Astrophysics},
+			publisher = {EDP Sciences},
+			author = {Neronov, A. and Savchenko, V. and Tramacere, A. and Meharga, M. and Ferrigno, C. and Paltani, S.},
+			year = {2021},
+			month = {jul},
+			pages = {A97}
+		}</citation>
+  </citations>
+</tool>
\ No newline at end of file