Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update tool MWA to 0.0.1+galaxy0 #81

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions tools/mwa/.shed.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
categories:
- Astronomy
description: MWA
homepage_url: null
long_description: MWA
name: mwa_astro_tool
owner: astroteam
remote_repository_url: https://github.com/esg-epfl-apc/tools-astro/tree/main/tools
type: unrestricted
91 changes: 91 additions & 0 deletions tools/mwa/Image.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
#!/usr/bin/env python
# coding: utf-8

# flake8: noqa

import json
import os
import shutil

import astropy.units as u
import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astroquery.skyview import SkyView
from oda_api.data_products import ImageDataProduct, PictureProduct
from oda_api.json import CustomJSONEncoder

src_name = "1ES 0229+200" # http://odahub.io/ontology#AstrophysicalObject
RA = 38.202562 # http://odahub.io/ontology#PointOfInterestRA
DEC = 20.288191 # http://odahub.io/ontology#PointOfInterestDEC
T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
Radius = 1.0 # http://odahub.io/ontology#AngleDegrees
pixsize = 0.01 # http://odahub.io/ontology#AngleDegrees
Frequency = "GLEAM 170-231 MHz" # http://odahub.io/ontology#String ; oda:allowed_value "GLEAM 72-103 MHz","GLEAM 103-134 MHz","GLEAM 139-170 MHz","GLEAM 170-231 MHz"

_galaxy_wd = os.getcwd()

with open("inputs.json", "r") as fd:
inp_dic = json.load(fd)
if "_data_product" in inp_dic.keys():
inp_pdic = inp_dic["_data_product"]
else:
inp_pdic = inp_dic

for vn, vv in inp_pdic.items():
if vn != "_selector":
globals()[vn] = type(globals()[vn])(vv)

pixels = int(2 * Radius / pixsize) + 1
Radius *= u.deg
pos = str(RA) + ", " + str(DEC)
pixels

hdul = SkyView.get_images(
position=pos, survey=[Frequency], pixels=pixels, radius=Radius
)

hdu = hdul[0]
hdu[0].header
wcs = WCS(hdu[0].header)

image = hdu[0].data

ax = plt.subplot(projection=wcs)
im = ax.imshow(image, origin="lower")
ax.coords.grid(True, color="white", ls="solid")
plt.colorbar(im, label="Jy/beam")
plt.savefig("Image.png", format="png", bbox_inches="tight")

hdu.writeto("Image.fits", overwrite=True)
bin_image = PictureProduct.from_file("Image.png")
fits_image = ImageDataProduct.from_fits_file("Image.fits")

picture = bin_image # http://odahub.io/ontology#ODAPictureProduct
image = fits_image # http://odahub.io/ontology#Image

# output gathering
_galaxy_meta_data = {}
_oda_outs = []
_oda_outs.append(("out_Image_picture", "picture_galaxy.output", picture))
_oda_outs.append(("out_Image_image", "image_galaxy.output", image))

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 ***")
117 changes: 117 additions & 0 deletions tools/mwa/Spectrum.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
#!/usr/bin/env python
# coding: utf-8

# flake8: noqa

import json
import os
import shutil

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.constants import h
from oda_api.data_products import ODAAstropyTable, PictureProduct
from oda_api.json import CustomJSONEncoder
from pyvo import registry # version >=1.4.1

src_name = "1ES 0229+200" # http://odahub.io/ontology#AstrophysicalObject
RA = 38.202562 # http://odahub.io/ontology#PointOfInterestRA
DEC = 20.288191 # http://odahub.io/ontology#PointOfInterestDEC
T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
Radius = 0.05 # http://odahub.io/ontology#AngleDegrees

_galaxy_wd = os.getcwd()

with open("inputs.json", "r") as fd:
inp_dic = json.load(fd)
if "_data_product" in inp_dic.keys():
inp_pdic = inp_dic["_data_product"]
else:
inp_pdic = inp_dic

for vn, vv in inp_pdic.items():
if vn != "_selector":
globals()[vn] = type(globals()[vn])(vv)

conesearch_radius = Radius # in degrees
conesearch_center = (RA, DEC)

# the catalogue name in VizieR
CATALOGUE = "VIII/100"
# each resource in the VO has an identifier, called ivoid. For vizier catalogs,
# the VO ids can be constructed like this:
catalogue_ivoid = f"ivo://CDS.VizieR/{CATALOGUE}"
# the actual query to the registry
voresource = registry.search(ivoid=catalogue_ivoid)[0]

conesearch_records = voresource.get_service("conesearch").search(
pos=conesearch_center,
sr=conesearch_radius,
)
conesearch_records

h_p = (h / u.s).to(u.eV).value # Planck constant in eV*s

conesearch_records.fieldnames
nu = []
F = []
F_err = []
for f in conesearch_records.fieldnames:
if (f[:4] == "Fint") and (f[4] != "w") and (f[4] != "f"):
nu.append(int(f[-3:]) * 1e6) # in Hz
F.append(conesearch_records[f][0] * nu[-1] * 1e-23)
F_err.append(conesearch_records["e_" + f][0] * nu[-1] * 1e-23)

E = h_p * np.array(nu)
plt.errorbar(E, F, F_err)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("$E$, eV")
plt.ylabel("$E F_E$, erg/cm$^2$s")
plt.savefig("Spectrum.png", format="png", bbox_inches="tight")

bin_image = PictureProduct.from_file("Spectrum.png")
from astropy.table import Table

data = [E, F, F_err]
names = ("E[eV]", "Flux[erg/cm2s]", "Flux_error[erg/cm2s]")
spec = ODAAstropyTable(Table(data, names=names))

picture_png = bin_image # http://odahub.io/ontology#ODAPictureProduct
spectrum_astropy_table = spec # http://odahub.io/ontology#ODAAstropyTable

# output gathering
_galaxy_meta_data = {}
_oda_outs = []
_oda_outs.append(
("out_Spectrum_picture_png", "picture_png_galaxy.output", picture_png)
)
_oda_outs.append(
(
"out_Spectrum_spectrum_astropy_table",
"spectrum_astropy_table_galaxy.output",
spectrum_astropy_table,
)
)

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 ***")
102 changes: 102 additions & 0 deletions tools/mwa/mwa_astro_tool.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
<tool id="mwa_astro_tool" name="MWA" version="0.0.1+galaxy0" profile="23.0">
<requirements>
<requirement type="package" version="8.28.0">ipython</requirement>
<requirement type="package" version="5.3.4">astropy</requirement>
<requirement type="package" version="1.2">gammapy</requirement>
<requirement type="package" version="3.9.2">matplotlib</requirement>
<requirement type="package" version="0.4.7">astroquery</requirement>
<requirement type="package" version="1.14.1">scipy</requirement>
<requirement type="package" version="1.2.28">oda-api</requirement>
<requirement type="package" version="7.16.4">nbconvert</requirement>
<requirement type="package" version="1.21.4">wget</requirement>
<requirement type="package" version="2.2.3">pandas</requirement>
<!--Requirements string 'nb2workflow[cwl,service,rdf,mmoda]>=1.3.30
' can't be converted automatically. Please add the galaxy/conda requirement manually or modify the requirements file!-->
</requirements>
<command detect_errors="exit_code">python '$__tool_directory__/${_data_product._selector}.py'</command>
<configfiles>
<inputs name="inputs" filename="inputs.json" data_style="paths" />
</configfiles>
<inputs>
<conditional name="_data_product">
<param name="_selector" type="select" label="Data Product">
<option value="Image" selected="true">Image</option>
<option value="Spectrum" selected="false">Spectrum</option>
</param>
<when value="Image">
<param name="src_name" type="text" value="1ES 0229+200" label="src_name" />
<param name="RA" type="float" value="38.202562" label="RA (unit: deg)" />
<param name="DEC" type="float" value="20.288191" label="DEC (unit: deg)" />
<param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" />
<param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" />
<param name="Radius" type="float" value="1.0" label="Radius (unit: deg)" />
<param name="pixsize" type="float" value="0.01" label="pixsize (unit: deg)" />
<param name="Frequency" type="select" label="Frequency">
<option value="GLEAM 103-134 MHz">GLEAM 103-134 MHz</option>
<option value="GLEAM 139-170 MHz">GLEAM 139-170 MHz</option>
<option value="GLEAM 170-231 MHz" selected="true">GLEAM 170-231 MHz</option>
<option value="GLEAM 72-103 MHz">GLEAM 72-103 MHz</option>
</param>
</when>
<when value="Spectrum">
<param name="src_name" type="text" value="1ES 0229+200" label="src_name" />
<param name="RA" type="float" value="38.202562" label="RA (unit: deg)" />
<param name="DEC" type="float" value="20.288191" label="DEC (unit: deg)" />
<param name="T1" type="text" value="2000-10-09T13:16:00.0" label="T1" />
<param name="T2" type="text" value="2022-10-10T13:16:00.0" label="T2" />
<param name="Radius" type="float" value="0.05" label="Radius (unit: deg)" />
</when>
</conditional>
</inputs>
<outputs>
<data label="${tool.name} -&gt; Image picture" name="out_Image_picture" format="auto" from_work_dir="picture_galaxy.output">
<filter>_data_product['_selector'] == 'Image'</filter>
</data>
<data label="${tool.name} -&gt; Image image" name="out_Image_image" format="auto" from_work_dir="image_galaxy.output">
<filter>_data_product['_selector'] == 'Image'</filter>
</data>
<data label="${tool.name} -&gt; Spectrum picture_png" name="out_Spectrum_picture_png" format="auto" from_work_dir="picture_png_galaxy.output">
<filter>_data_product['_selector'] == 'Spectrum'</filter>
</data>
<data label="${tool.name} -&gt; Spectrum spectrum_astropy_table" name="out_Spectrum_spectrum_astropy_table" format="auto" from_work_dir="spectrum_astropy_table_galaxy.output">
<filter>_data_product['_selector'] == 'Spectrum'</filter>
</data>
</outputs>
<tests>
<test expect_num_outputs="2">
<conditional name="_data_product">
<param name="_selector" value="Image" />
<param name="src_name" value="1ES 0229+200" />
<param name="RA" value="38.202562" />
<param name="DEC" value="20.288191" />
<param name="T1" value="2000-10-09T13:16:00.0" />
<param name="T2" value="2022-10-10T13:16:00.0" />
<param name="Radius" value="1.0" />
<param name="pixsize" value="0.01" />
<param name="Frequency" value="GLEAM 170-231 MHz" />
</conditional>
<assert_stdout>
<has_text text="*** Job finished successfully ***" />
</assert_stdout>
</test>
<test expect_num_outputs="2">
<conditional name="_data_product">
<param name="_selector" value="Spectrum" />
<param name="src_name" value="1ES 0229+200" />
<param name="RA" value="38.202562" />
<param name="DEC" value="20.288191" />
<param name="T1" value="2000-10-09T13:16:00.0" />
<param name="T2" value="2022-10-10T13:16:00.0" />
<param name="Radius" value="0.05" />
</conditional>
<assert_stdout>
<has_text text="*** Job finished successfully ***" />
</assert_stdout>
</test>
</tests>
<help>This is help for MWA
</help>
<citations>
<citation type="doi">10.1093/mnras/stw2337</citation>
</citations>
</tool>
Loading