Mercurial > repos > astroteam > grb_detection_astro_tool
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} -> detectgrb lc" name="out_detectgrb_lc" format="auto" from_work_dir="lc_galaxy.output" /> + <data label="${tool.name} -> detectgrb detection_comment" name="out_detectgrb_detection_comment" format="auto" from_work_dir="detection_comment_galaxy.output" /> + <data label="${tool.name} -> detectgrb image_output" name="out_detectgrb_image_output" format="auto" from_work_dir="image_output_galaxy.output" /> + <data label="${tool.name} -> 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; 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