Mercurial > repos > astroteam > plot_tools_astro_tool

comparison light_curve.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 "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 ***")