Mercurial > repos > astroteam > plot_tools_astro_tool

comparison spectrum.py @ 0:2b1759ccaa8b draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 ***")