Mercurial > repos > astroteam > plot_tools_astro_tool
comparison light_curve.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 "comma", "tab", "space", "whitespace", "semicolon" | |
19 column = "T" # 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 use_quantile_values = False # http://odahub.io/ontology#Boolean | |
25 nbins = 15 # http://odahub.io/ontology#Integer | |
26 xlabel = "time, s" # http://odahub.io/ontology#String | |
27 ylabel = "Ncounts" # http://odahub.io/ontology#String | |
28 plot_mode = "flux" # http://odahub.io/ontology#String ; oda:allowed_value "counts", "flux" | |
29 | |
30 _galaxy_wd = os.getcwd() | |
31 | |
32 with open("inputs.json", "r") as fd: | |
33 inp_dic = json.load(fd) | |
34 if "C_data_product_" in inp_dic.keys(): | |
35 inp_pdic = inp_dic["C_data_product_"] | |
36 else: | |
37 inp_pdic = inp_dic | |
38 fn = str(inp_pdic["fn"]) | |
39 skiprows = int(inp_pdic["skiprows"]) | |
40 sep = str(inp_pdic["sep"]) | |
41 column = str(inp_pdic["column"]) | |
42 weight_col = str(inp_pdic["weight_col"]) | |
43 binning = str(inp_pdic["binning"]) | |
44 minval = float(inp_pdic["minval"]) | |
45 maxval = float(inp_pdic["maxval"]) | |
46 use_quantile_values = bool(inp_pdic["use_quantile_values"]) | |
47 nbins = int(inp_pdic["nbins"]) | |
48 xlabel = str(inp_pdic["xlabel"]) | |
49 ylabel = str(inp_pdic["ylabel"]) | |
50 plot_mode = str(inp_pdic["plot_mode"]) | |
51 | |
52 import matplotlib.pyplot as plt | |
53 import numpy as np | |
54 import pandas as pd | |
55 | |
56 if plot_mode != "counts" and ylabel == "Ncounts": | |
57 ylabel = plot_mode # replace default value | |
58 | |
59 assert minval >= 0 or not use_quantile_values | |
60 assert maxval >= 0 or not use_quantile_values | |
61 assert minval <= 1 or not use_quantile_values | |
62 assert maxval <= 1 or not use_quantile_values | |
63 assert minval < maxval or minval == 0 or maxval == 0 | |
64 | |
65 separators = { | |
66 "tab": "\t", | |
67 "comma": ",", | |
68 "semicolon": ";", | |
69 "whitespace": "\s+", | |
70 "space": " ", | |
71 } | |
72 | |
73 df = None | |
74 | |
75 if sep == "auto": | |
76 for name, s in separators.items(): | |
77 try: | |
78 df = pd.read_csv(fn, sep=s, index_col=False, skiprows=skiprows) | |
79 if len(df.columns) > 2: | |
80 sep = s | |
81 print("Detected separator: ", name) | |
82 break | |
83 except Exception as e: | |
84 print("Separator ", s, " failed", e) | |
85 assert sep != "auto", "Failed to find valid separator" | |
86 | |
87 if df is None: | |
88 df = pd.read_csv(fn, sep=separators[sep], index_col=False) | |
89 | |
90 df.columns | |
91 | |
92 def weighted_quantile( | |
93 values, quantiles, sample_weight=None, values_sorted=False, old_style=False | |
94 ): | |
95 """Very close to numpy.percentile, but supports weights. | |
96 NOTE: quantiles should be in [0, 1]! | |
97 :param values: numpy.array with data | |
98 :param quantiles: array-like with many quantiles needed | |
99 :param sample_weight: array-like of the same length as `array` | |
100 :param values_sorted: bool, if True, then will avoid sorting of initial array | |
101 :param old_style: if True, will correct output to be consistent with numpy.percentile. | |
102 :return: numpy.array with computed quantiles. | |
103 """ | |
104 values = np.array(values) | |
105 quantiles = np.array(quantiles) | |
106 if sample_weight is None: | |
107 sample_weight = np.ones(len(values)) | |
108 sample_weight = np.array(sample_weight) | |
109 assert np.all(quantiles >= 0) and np.all( | |
110 quantiles <= 1 | |
111 ), "quantiles should be in [0, 1]" | |
112 | |
113 if not values_sorted: | |
114 sorter = np.argsort(values) | |
115 values = values[sorter] | |
116 sample_weight = sample_weight[sorter] | |
117 | |
118 weighted_quantiles = np.cumsum(sample_weight) - 0.5 * sample_weight | |
119 if old_style: | |
120 # To be convenient with np.percentile | |
121 weighted_quantiles -= weighted_quantiles[0] | |
122 weighted_quantiles /= weighted_quantiles[-1] | |
123 else: | |
124 weighted_quantiles /= np.sum(sample_weight) | |
125 return np.interp(quantiles, weighted_quantiles, values) | |
126 | |
127 def read_data(df, colname, optional=False): | |
128 for i, c in enumerate(df.columns): | |
129 if colname == f"c{i+1}": | |
130 print(colname, c) | |
131 return df[c].values | |
132 elif colname == c: | |
133 print(colname, c) | |
134 return df[c].values | |
135 | |
136 assert optional, colname + " column not found" | |
137 return None | |
138 | |
139 delays = read_data(df, column) | |
140 weights = read_data(df, weight_col, optional=True) | |
141 if weights is None: | |
142 weights = np.ones_like(delays) | |
143 | |
144 if binning != "linear": | |
145 min_positive_val = np.min(delays[delays > 0]) | |
146 delays[delays <= 0] = ( | |
147 min_positive_val # replace zero delays with minimal positive value | |
148 ) | |
149 | |
150 if use_quantile_values: | |
151 minval, maxval = weighted_quantile( | |
152 delays, [minval, maxval], sample_weight=weights | |
153 ) | |
154 if minval == maxval: | |
155 print("ignoreing minval and maxval (empty range)") | |
156 minval = np.min(delays) | |
157 maxval = np.max(delays) | |
158 else: | |
159 if minval == 0: | |
160 minval = np.min(delays) | |
161 if maxval == 0: | |
162 maxval = np.max(delays) | |
163 | |
164 if minval == maxval: | |
165 print("correcting minval and maxval (empty range)") | |
166 maxval = minval * 1.1 if minval > 0 else 1e-100 | |
167 | |
168 from numpy import log10 | |
169 | |
170 if binning == "linear": | |
171 bins = np.linspace(minval, maxval, nbins + 1) | |
172 else: | |
173 bins = np.logspace(log10(minval), log10(maxval), nbins + 1) | |
174 bins | |
175 | |
176 if plot_mode == "flux" and binning == "logarithmic": | |
177 weights = weights / delays | |
178 | |
179 plt.figure() | |
180 h = plt.hist(delays, weights=weights, bins=bins) | |
181 | |
182 if binning == "logarithmic": | |
183 plt.xscale("log") | |
184 plt.yscale("log") | |
185 plt.xlabel(xlabel) | |
186 plt.ylabel(ylabel) | |
187 plt.savefig("Histogram.png", format="png", dpi=150) | |
188 hist_counts = h[0] | |
189 hist_bins = h[1] | |
190 hist_mins = hist_bins[:-1] | |
191 hist_maxs = hist_bins[1:] | |
192 | |
193 from astropy.table import Table | |
194 from oda_api.data_products import ODAAstropyTable, PictureProduct | |
195 | |
196 names = ("bins_min", "bins_max", "counts") | |
197 res = ODAAstropyTable(Table([hist_mins, hist_maxs, hist_counts], names=names)) | |
198 | |
199 plot = PictureProduct.from_file("Histogram.png") | |
200 | |
201 histogram_data = res # http://odahub.io/ontology#ODAAstropyTable | |
202 histogram_picture = plot # http://odahub.io/ontology#ODAPictureProduct | |
203 | |
204 # output gathering | |
205 _galaxy_meta_data = {} | |
206 _oda_outs = [] | |
207 _oda_outs.append( | |
208 ( | |
209 "out_light_curve_histogram_data", | |
210 "histogram_data_galaxy.output", | |
211 histogram_data, | |
212 ) | |
213 ) | |
214 _oda_outs.append( | |
215 ( | |
216 "out_light_curve_histogram_picture", | |
217 "histogram_picture_galaxy.output", | |
218 histogram_picture, | |
219 ) | |
220 ) | |
221 | |
222 for _outn, _outfn, _outv in _oda_outs: | |
223 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) | |
224 if isinstance(_outv, str) and os.path.isfile(_outv): | |
225 shutil.move(_outv, _galaxy_outfile_name) | |
226 _galaxy_meta_data[_outn] = {"ext": "_sniff_"} | |
227 elif getattr(_outv, "write_fits_file", None): | |
228 _outv.write_fits_file(_galaxy_outfile_name) | |
229 _galaxy_meta_data[_outn] = {"ext": "fits"} | |
230 elif getattr(_outv, "write_file", None): | |
231 _outv.write_file(_galaxy_outfile_name) | |
232 _galaxy_meta_data[_outn] = {"ext": "_sniff_"} | |
233 else: | |
234 with open(_galaxy_outfile_name, "w") as fd: | |
235 json.dump(_outv, fd, cls=CustomJSONEncoder) | |
236 _galaxy_meta_data[_outn] = {"ext": "json"} | |
237 | |
238 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd: | |
239 json.dump(_galaxy_meta_data, fd) | |
240 print("*** Job finished successfully ***") |