Mercurial > repos > astroteam > plot_tools_astro_tool
comparison spectrum.py @ 0:2b1759ccaa8b draft default tip
planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit f28a8cb73a7f3053eac92166867a48b3d4af28fd
author | astroteam |
---|---|
date | Fri, 25 Apr 2025 21:48:27 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:2b1759ccaa8b |
---|---|
1 #!/usr/bin/env python | |
2 # coding: utf-8 | |
3 | |
4 #!/usr/bin/env python | |
5 | |
6 # This script is generated with nb2galaxy | |
7 | |
8 # flake8: noqa | |
9 | |
10 import json | |
11 import os | |
12 import shutil | |
13 | |
14 from oda_api.json import CustomJSONEncoder | |
15 | |
16 fn = "data.tsv" # oda:POSIXPath | |
17 skiprows = 0 # http://odahub.io/ontology#Integer | |
18 sep = "whitespace" # http://odahub.io/ontology#String ; oda:allowed_value "auto", "comma", "tab", "whitespace", "semicolon" | |
19 column = "c1" # http://odahub.io/ontology#String | |
20 weight_col = "" # http://odahub.io/ontology#String | |
21 binning = "logarithmic" # http://odahub.io/ontology#String ; oda:allowed_value "linear","logarithmic" | |
22 minval = 0 # http://odahub.io/ontology#Float | |
23 maxval = 0 # http://odahub.io/ontology#Float | |
24 nbins = 15 # http://odahub.io/ontology#Integer | |
25 xlabel = "Energy, [eV]" # http://odahub.io/ontology#String | |
26 ylabel = "Flux E^2, [eV]" # http://odahub.io/ontology#String | |
27 spec_power = 2.0 # http://odahub.io/ontology#Float | |
28 | |
29 _galaxy_wd = os.getcwd() | |
30 | |
31 with open("inputs.json", "r") as fd: | |
32 inp_dic = json.load(fd) | |
33 if "C_data_product_" in inp_dic.keys(): | |
34 inp_pdic = inp_dic["C_data_product_"] | |
35 else: | |
36 inp_pdic = inp_dic | |
37 fn = str(inp_pdic["fn"]) | |
38 skiprows = int(inp_pdic["skiprows"]) | |
39 sep = str(inp_pdic["sep"]) | |
40 column = str(inp_pdic["column"]) | |
41 weight_col = str(inp_pdic["weight_col"]) | |
42 binning = str(inp_pdic["binning"]) | |
43 minval = float(inp_pdic["minval"]) | |
44 maxval = float(inp_pdic["maxval"]) | |
45 nbins = int(inp_pdic["nbins"]) | |
46 xlabel = str(inp_pdic["xlabel"]) | |
47 ylabel = str(inp_pdic["ylabel"]) | |
48 spec_power = float(inp_pdic["spec_power"]) | |
49 | |
50 import matplotlib.pyplot as plt | |
51 import numpy as np | |
52 import pandas as pd | |
53 | |
54 separators = { | |
55 "tab": "\t", | |
56 "comma": ",", | |
57 "semicolon": ";", | |
58 "whitespace": "\s+", | |
59 "space": " ", | |
60 } | |
61 | |
62 df = None | |
63 | |
64 if sep == "auto": | |
65 for name, s in separators.items(): | |
66 try: | |
67 df = pd.read_csv(fn, sep=s, index_col=False, skiprows=skiprows) | |
68 if len(df.columns) > 2: | |
69 sep = s | |
70 print("Detected separator: ", name) | |
71 break | |
72 except Exception as e: | |
73 print("Separator ", s, " failed", e) | |
74 assert sep != "auto", "Failed to find valid separator" | |
75 | |
76 if df is None: | |
77 df = pd.read_csv(fn, sep=separators[sep], index_col=False) | |
78 | |
79 df.columns | |
80 | |
81 def read_data(df, colname, optional=False): | |
82 for i, c in enumerate(df.columns): | |
83 if colname == f"c{i+1}": | |
84 print(colname, c) | |
85 return df[c].values | |
86 elif colname == c: | |
87 print(colname, c) | |
88 return df[c].values | |
89 | |
90 assert optional, colname + " column not found" | |
91 return None | |
92 | |
93 values = read_data(df, column) | |
94 weights = read_data(df, weight_col, optional=True) | |
95 if weights is None: | |
96 weights = np.ones_like(values) | |
97 | |
98 values, weights | |
99 | |
100 from numpy import log10 | |
101 | |
102 if minval == 0: | |
103 minval = np.min(values) | |
104 | |
105 if maxval == 0: | |
106 maxval = np.max(values) | |
107 | |
108 if binning == "linear": | |
109 bins = np.linspace(minval, maxval, nbins + 1) | |
110 else: | |
111 bins = np.logspace(log10(minval), log10(maxval), nbins + 1) | |
112 bins | |
113 | |
114 bin_val, _ = np.histogram(values, weights=weights, bins=bins) | |
115 len(bin_val), len(bins) | |
116 bin_width = bins[1:] - bins[:-1] | |
117 flux = bin_val / bin_width | |
118 if binning == "linear": | |
119 spec_point = 0.5 * (bins[1:] + bins[:-1]) | |
120 else: | |
121 spec_point = np.sqrt(bins[1:] * bins[:-1]) | |
122 | |
123 plt.figure() | |
124 h = plt.plot(spec_point, flux * spec_point**spec_power) | |
125 | |
126 if binning == "logarithmic": | |
127 plt.xscale("log") | |
128 plt.yscale("log") | |
129 | |
130 plt.xlabel(xlabel) | |
131 plt.ylabel(ylabel) | |
132 plt.savefig("spectrum.png", format="png", dpi=150) | |
133 | |
134 from astropy.table import Table | |
135 from oda_api.data_products import ODAAstropyTable, PictureProduct | |
136 | |
137 names = ("bins_min", "bins_max", "flux") | |
138 | |
139 res = ODAAstropyTable(Table([bins[:-1], bins[1:], flux], names=names)) | |
140 | |
141 plot = PictureProduct.from_file("spectrum.png") | |
142 | |
143 histogram_data = res # http://odahub.io/ontology#ODAAstropyTable | |
144 histogram_picture = plot # http://odahub.io/ontology#ODAPictureProduct | |
145 | |
146 # output gathering | |
147 _galaxy_meta_data = {} | |
148 _oda_outs = [] | |
149 _oda_outs.append( | |
150 ( | |
151 "out_spectrum_histogram_data", | |
152 "histogram_data_galaxy.output", | |
153 histogram_data, | |
154 ) | |
155 ) | |
156 _oda_outs.append( | |
157 ( | |
158 "out_spectrum_histogram_picture", | |
159 "histogram_picture_galaxy.output", | |
160 histogram_picture, | |
161 ) | |
162 ) | |
163 | |
164 for _outn, _outfn, _outv in _oda_outs: | |
165 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) | |
166 if isinstance(_outv, str) and os.path.isfile(_outv): | |
167 shutil.move(_outv, _galaxy_outfile_name) | |
168 _galaxy_meta_data[_outn] = {"ext": "_sniff_"} | |
169 elif getattr(_outv, "write_fits_file", None): | |
170 _outv.write_fits_file(_galaxy_outfile_name) | |
171 _galaxy_meta_data[_outn] = {"ext": "fits"} | |
172 elif getattr(_outv, "write_file", None): | |
173 _outv.write_file(_galaxy_outfile_name) | |
174 _galaxy_meta_data[_outn] = {"ext": "_sniff_"} | |
175 else: | |
176 with open(_galaxy_outfile_name, "w") as fd: | |
177 json.dump(_outv, fd, cls=CustomJSONEncoder) | |
178 _galaxy_meta_data[_outn] = {"ext": "json"} | |
179 | |
180 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd: | |
181 json.dump(_galaxy_meta_data, fd) | |
182 print("*** Job finished successfully ***") |