diff --git a/README.md b/README.md index 651ee01..df4e55f 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # uds -Содержит код для обработки данных поля UDS. +Содержит код для отрисовки данных монитора ART-XC. # Установка @@ -10,26 +10,23 @@ mkdir ~/work; cd ~/work python -m venv venv source ./venv/bin/activate.csh ``` -или запускаем другое виртуальное окружение, например: - -``` conda activate ciao-4.15 ``` ## Шаг 2. Клонируем проект из репозитория ``` mkdir ~/work; cd ~/work -git clone git@danko.iki.rssi.ru:erosita/uds.git +git clone git@heagit.cosmos.ru:artxc/monitor.git ``` ## Шаг 3. Устанавливаем проект в ваше виртуальное окружение ``` -cd uds -pip install --editable uds/ +cd monitor +pip install --editable monitor/ ``` Обратите внимание на параметр **--editable**, он позволяет вам редактировать исходный код данного пакета и сразу его выполнять. Если вы не планируете модифицировать локальную копию кода, можете убрать этот параметр. После работы можете удалить проект: -``` pip uninstall uds ``` +``` pip uninstall monitor ``` ## Работа с данными diff --git a/monitor/monitor/config.py b/monitor/monitor/config.py index e69de29..7f6db14 100644 --- a/monitor/monitor/config.py +++ b/monitor/monitor/config.py @@ -0,0 +1,3 @@ +threshold=1.1 + +download_link='https://monitor.srg.cosmos.ru/download?type=orig&range=60.0-120.0' diff --git a/monitor/monitor/utils.py b/monitor/monitor/utils.py index d455957..70ee51f 100644 --- a/monitor/monitor/utils.py +++ b/monitor/monitor/utils.py @@ -1,6 +1,17 @@ import os import sys import csv +import numpy as np + +import datetime +import dateutil +from scipy.optimize import curve_fit +from scipy.ndimage import gaussian_filter1d +from scipy import signal +from scipy.stats import median_abs_deviation +import astropy.io.fits as pf +from astropy.time import Time +from datetime import datetime, timedelta def create_folder(folder): if not (os.path.exists(folder)): @@ -10,3 +21,153 @@ def remove_file(filename): if(os.path.isfile(filename)==True): os.remove(filename) +def polynomial_func(x, *coefficients): + """General polynomial function""" + return sum(coef * x**i for i, coef in enumerate(coefficients)) + +def detect_outbursts_threshold(times, fluxes, threshold=2.0, min_duration=3): + """ + Detect outbursts using a threshold-based approach. + + Parameters: + times: array of time values + fluxes: array of flux values + threshold: number of standard deviations above median to consider as outburst + min_duration: minimum number of consecutive points to consider as outburst + + Returns: + outburst_intervals: list of (start_time, end_time) tuples + """ + # Calculate median and standard deviation + median_flux = np.median(fluxes) + mad = median_abs_deviation(fluxes) + print("median_flux",median_flux,"mad",mad) + + # Create a threshold value + threshold_value = median_flux + threshold * mad + + # Identify points above threshold + above_threshold = fluxes > threshold_value + + # Find consecutive points above threshold + outburst_indices = [] + in_outburst = False + start_idx = None + + for i, is_above in enumerate(above_threshold): + if is_above and not in_outburst: + in_outburst = True + start_idx = i + elif not is_above and in_outburst: + in_outburst = False + # Check if the outburst duration meets the minimum requirement + if i - start_idx >= min_duration: + outburst_indices.append((start_idx, i-1)) + start_idx = None + + # Handle case where outburst continues until end of data + if in_outburst and len(fluxes) - start_idx >= min_duration: + outburst_indices.append((start_idx, len(fluxes)-1)) + + # Convert indices to time intervals + outburst_intervals = [(times[start], times[end]) for start, end in outburst_indices] + + return outburst_intervals, median_flux,mad + +def detect_outbursts(times, fluxes, threshold=2.0, min_duration=3): + """ + Detect outbursts using a threshold-based approach. + + Parameters: + times: array of time values + fluxes: array of flux values + threshold: number of standard deviations above median to consider as outburst + min_duration: minimum number of consecutive points to consider as outburst + + Returns: + outburst_intervals: list of (start_time, end_time) tuples + """ + + # Identify points above threshold + above_threshold = fluxes > threshold + + # Find consecutive points above threshold + outburst_indices = [] + in_outburst = False + start_idx = None + peak_idx = None + peak_flux = 0 + + for i, is_above in enumerate(above_threshold): + if is_above and not in_outburst: + in_outburst = True + start_idx = i + elif not is_above and in_outburst: + # fist point after outburst + in_outburst = False + # Check if the outburst duration meets the minimum requirement + if i - start_idx >= min_duration: + outburst_indices.append((start_idx, i-1, peak_idx)) + start_idx = None + peak_idx = None + peak_flux = 0 + if is_above and in_outburst: + if peak_flux < fluxes[i]: + peak_flux = fluxes[i] + peak_idx = i + + + # Handle case where outburst continues until end of data + if in_outburst and len(fluxes) - start_idx >= min_duration: + outburst_indices.append((start_idx, len(fluxes)-1, peak_idx)) # not tested + + + # Convert indices to time intervals + outburst_intervals = [(times[start], times[end], times[peak], fluxes[peak]) for start, end, peak in outburst_indices] + + + return outburst_intervals + +def get_decimal_year(dt): + day_of_year = dt.timetuple().tm_yday+dt.timetuple().tm_hour/24+dt.timetuple().tm_min/60/24+dt.timetuple().tm_sec/60/60/24 # returns 1 for January 1st + is_leap = (dt.year % 4 == 0 and dt.year % 100 != 0) or (dt.year % 400 == 0) + days_in_year = 366 if is_leap else 365 + return dt.year+(day_of_year/days_in_year) + +# AI generated: +def decimal_year_to_date(decimal_year): + """ + Converts a decimal year (e.g., 2023.5) to a datetime object. + Accounts for leap years. + """ + year = int(decimal_year) + fractional_part = decimal_year - year + + # Determine if it's a leap year to get the correct number of days + is_leap = (year % 4 == 0 and year % 100 != 0) or (year % 400 == 0) + days_in_year = 366 if is_leap else 365 + + # Calculate the number of days into the year + days_into_year = fractional_part * days_in_year + + # Create a datetime object for January 1st of the given year + start_of_year = datetime(year, 1, 1) + + # Add the calculated number of days to get the final date + # Subtract 1 day because the fractional part represents days *after* Jan 1st + result_date = start_of_year + timedelta(days=days_into_year) + + return result_date + +def save_tex_table(outbursts, output=None): + fmt="%d %b %H:%M" + with open(output, 'w') as file: + for start, end, peak, peak_flux in outbursts: + date_start = decimal_year_to_date(start) + date_end = decimal_year_to_date(end) + date_peak = decimal_year_to_date(peak) + print_year = date_start.strftime("%Y") + print_start = date_start.strftime(fmt) + print_end = date_end.strftime(fmt) + print_peak = date_peak.strftime(fmt) + file.write(f'{print_year} & {print_start} & {print_peak} & {peak_flux:.2f} & {print_end} \\\\\n') diff --git a/scripts/download.py b/scripts/download.py new file mode 100755 index 0000000..daf6fcc --- /dev/null +++ b/scripts/download.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + plot.py + + +НАЗНАЧЕНИЕ: + + Простой скрипт для отрисовки файла, выгруженного с СРГ L2 монитора https://monitor.srg.cosmos.ru/ + +ВЫЗОВ: + + conda activate + ./plot.py + + +УПРАВЛЕНИЕ: + + Название файла надо вставить внутрь скрипта (ищите default.csv) + +ПАРАМЕТРЫ: + + N/A + + +ВЫВОД: + + Файл monitor.png записывается в текущую директорию + + +ИСТОРИЯ: + Роман Кривонос, ИКИ РАН, krivonos@cosmos.ru + Декабрь 2024 + +""" +import urllib +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import ticker +import pandas as pd +import datetime +import dateutil +from os.path import dirname +import inspect + +import requests + +import monitor +from monitor.config import * +from monitor.utils import * + +""" find root dir """ +root_path=dirname(dirname(dirname(inspect.getfile(monitor)))) +print("Monitor root path: {}".format(root_path)) + + +products_dir=root_path+'/products/' +data_dir=root_path+'/data/' + +create_folder(products_dir) + +fn='default.csv' + +local_filename=data_dir+fn + +url = urllib.parse.quote(download_link, safe='/', encoding=None, errors=None) + +response = requests.get(url) +print(response.status_code) +if response.status_code == 200: + # Open the local file in binary write mode ('wb') and save the content + with open(local_filename, 'wb') as file: + file.write(response.content) + print(f"File '{local_filename}' downloaded successfully.") +else: + print(f"Failed to download file. Status code: {response.status_code}") diff --git a/scripts/gti.py b/scripts/gti.py new file mode 100755 index 0000000..c6a6de1 --- /dev/null +++ b/scripts/gti.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + gti.py + + +НАЗНАЧЕНИЕ: + + Cкрипт для создания интервалов хорошего времени по данным из файла, сбинированного по 60 минут программой create_binned_filed.py (на основе оригинального файла lc_orig_60-120.csv, выгруженного с СРГ L2 монитора https://monitor.srg.cosmos.ru/ ) + Скрипт аппроксимирует кривую блеска полиномом, чтобы убрать синусоидальный тренд, ищет времена вспышек и возвращает интервалы времени без вспышек в офрмате секунд с начала 1 января 2000 года (Москва) (файл gti_sec.csv) и в формате дестичного года (файл gti_dyear.csv). + +ВЫЗОВ: + + conda activate + ./gti.py + + +УПРАВЛЕНИЕ: + + Название файла надо вставить внутрь скрипта (ищите lc_bin_60min_60-120.csv) + +ПАРАМЕТРЫ: + + N/A + + +ВЫВОД: + Файл с графиком monitor.png записывается в текущую директорию. + Файл с интервалами времени gti.csv записывается в текущую директорию. + + +ИСТОРИЯ: + Филиппова Екатерина, ИКИ РАН, kate@cosmos.ru + Сентябрь, 2025 + +""" + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import ticker +import pandas as pd +import datetime +import dateutil +from os.path import dirname +import inspect + +from scipy.optimize import curve_fit +from scipy.ndimage import gaussian_filter1d +from scipy import signal +from scipy.stats import median_abs_deviation +import astropy.io.fits as pf +from astropy.time import Time + + +import monitor +from monitor.config import * +from monitor.utils import * + +""" find root dir """ +root_path=dirname(dirname(dirname(inspect.getfile(monitor)))) +print("Monitor root path: {}".format(root_path)) + + +products_dir=root_path+'/products/' +data_dir=root_path+'/data/' + +create_folder(products_dir) + +fig, ax = plt.subplots(figsize=(9, 5), dpi=100) +ax.grid(visible=True,linestyle='dotted', ) +ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12)) +ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2)) +ax.tick_params(axis="both", width=1, labelsize=14) + +#plt.xlim([2020,2025]) +#plt.ylim([2,12]) +plt.ylabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal') +plt.xlabel('Year',fontsize=14, fontweight='normal') + +for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + +cl='black' + +input_file='default.csv' + +# Сюда надо вставить название файла, скачанного с https://monitor.srg.cosmos.ru/ +df = pd.read_csv(data_dir+input_file) + +tm=[] +rate=[] +error=[] +for index, row in df.iterrows(): + dt = dateutil.parser.parse(row['timestamp']) + dyear = get_decimal_year(dt) + tm.append(dyear) + rate.append(float(row['value_60.0-120.0'])) + error.append(float(row['error_60.0-120.0'])) + +rate=np.array(rate) +tm=np.array(tm) +tmp=np.where(rate>0,1,0) +rate=np.compress(tmp,rate) +tm=np.compress(tmp,tm) +error=np.compress(tmp,error) + +tmp=np.where(rate < 4.3,1,0) +tm_fit=np.compress(tmp,tm) +rate_fit=np.compress(tmp,rate) + +# fit polinom of degree +degree = 6 +#initial_guess = [1] * (degree + 1) # Initial guess for coefficients +#initial_guess = [-3.67265411e+08, 5.44675290e+05, -2.69261063e+02, 4.43698213e-02] +initial_guess = [ -6.10911582e+07, 3.02005656e+04, 1.49297834e+01, 6.53394516e-09, -3.64859359e-06, -1.80368300e-09, 8.91651441e-13] + +params, covariance = curve_fit(polynomial_func, tm_fit, rate_fit, p0=initial_guess) +print("curve_fit:",params) + +# Generate fitted curve (original bin) +y_fit = polynomial_func(tm, *params) + +# Delaem interpolyaciyu v 100 raz +year_min = min(tm) +year_max = max(tm) +year_step = (year_max - year_min)/len(tm) +year_sim = np.arange(year_min, year_max, year_step/100) +rate_interp = np.interp(year_sim, tm, rate) + +# Generate fitted curve (interpolated) +y_fit_interp = polynomial_func(year_sim, *params) + +plt.plot(tm, rate, color=cl, linewidth=1, linestyle='solid') +plt.plot(tm, y_fit, color="red", linewidth=2, linestyle='solid', label=f'Degree {degree} fit') +plt.ylabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal') +plt.xlabel('Year',fontsize=14, fontweight='normal') +fig.savefig(products_dir+"bestfit.png", bbox_inches='tight') +plt.show() +fig.clear() # Clears the specific figure object + +fig, ax = plt.subplots(figsize=(9, 5), dpi=100) +ax.grid(visible=True,linestyle='dotted', ) +ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12)) +ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2)) +ax.tick_params(axis="both", width=1, labelsize=14) + +#plt.xlim([2020,2025]) +#plt.ylim([0,200]) + +for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + + +resid=rate/y_fit +resid_interp=rate_interp/y_fit_interp +#smoothed_resid = gaussian_filter1d(resid, sigma=2) +outbursts = detect_outbursts(year_sim, resid_interp, threshold=threshold, min_duration=1) + + +plt.plot(tm, resid, color="black", linewidth=1, linestyle='solid',label="Residuals") +#plt.plot(tm, smoothed_resid, color="red", linewidth=1, linestyle='solid',label="Residuals") +plt.axhline(y = threshold, color = 'b', linewidth=1 , linestyle='dashed', label = 'Threshold') +plt.ylabel('Ratio',fontsize=14, fontweight='normal') +plt.xlabel('Year',fontsize=14, fontweight='normal') +fig.savefig(products_dir+"resid.png", bbox_inches='tight') + + +# Testing +""" +for start,end,peak,peak_flux in outbursts: + print(start,end,peak,peak_flux) + plt.axvline(x = start, color = 'r', linewidth=1 , linestyle='solid') + plt.axvline(x = end, color = 'g', linewidth=1 , linestyle='solid') + plt.axvline(x = peak, color = 'b', linewidth=1 , linestyle='solid') +""" +plt.show() + + +save_tex_table(outbursts,output=products_dir+'flares.tex') + +tstart=[tm[0]] +tstop=[] +for start,end,peak,peak_flux in outbursts: + if end != tm[-1]: + tstop.append(start) + tstart.append(end) + else: + tstop.append(start) + +if outbursts[-1][1] != tm[-1]: + tstop.append(tm[-1]) + +df = pd.DataFrame({'tstart': tstart,'tstop': tstop}) +df.to_csv(products_dir+'gti_dyear.csv', index=False) + +tstart_obj = Time(tstart, format='decimalyear') +tstop_obj = Time(tstop,format='decimalyear') +tstart_sec=(tstart_obj.mjd-51543.75)*86400 +tstop_sec=(tstop_obj.mjd-51543.75)*86400 + +df = pd.DataFrame({'tstart': tstart_sec,'tstop': tstop_sec}) +df.to_csv(products_dir+'gti_sec.csv', index=False) + + diff --git a/scripts/plot.py b/scripts/plot.py new file mode 100755 index 0000000..a7e45d2 --- /dev/null +++ b/scripts/plot.py @@ -0,0 +1,111 @@ +#!/usr/bin/env python +""" +НАЗВАНИЕ: + + plot.py + + +НАЗНАЧЕНИЕ: + + Простой скрипт для отрисовки файла, выгруженного с СРГ L2 монитора https://monitor.srg.cosmos.ru/ + +ВЫЗОВ: + + conda activate + ./plot.py + + +УПРАВЛЕНИЕ: + + Название файла надо вставить внутрь скрипта (ищите default.csv) + +ПАРАМЕТРЫ: + + N/A + + +ВЫВОД: + + Файл monitor.png записывается в текущую директорию + + +ИСТОРИЯ: + Роман Кривонос, ИКИ РАН, krivonos@cosmos.ru + Декабрь 2024 + +""" + +import matplotlib.pyplot as plt +import numpy as np +from matplotlib import ticker +import pandas as pd +import datetime +import dateutil +from os.path import dirname +import inspect + +import monitor +from monitor.config import * +from monitor.utils import * + +""" find root dir """ +root_path=dirname(dirname(dirname(inspect.getfile(monitor)))) +print("Monitor root path: {}".format(root_path)) + + +products_dir=root_path+'/products/' +data_dir=root_path+'/data/' + +create_folder(products_dir) + +fig, ax = plt.subplots(figsize=(9, 5), dpi=100) + + +#ax.set(xlabel='time (s)', ylabel='voltage (mV)', +# title='About as simple as it gets, folks') + +ax.grid(visible=True,linestyle='dotted', ) +ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12)) +ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2)) +ax.tick_params(axis="both", width=1, labelsize=14) + +#ax.grid() +#plt.xlim([2020,2025]) +#plt.ylim([2,12]) +ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12)) +ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2)) +ax.tick_params(axis="both", width=1, labelsize=14) + +for axis in ['top','bottom','left','right']: + ax.spines[axis].set_linewidth(1) + +cl='black' + +input_file='default.csv' + +# Сюда надо вставить название файла, скачанного с https://monitor.srg.cosmos.ru/ +df = pd.read_csv(data_dir+input_file) + + +tm=[] +rate=[] + +for index, row in df.iterrows(): + dt = dateutil.parser.parse(row['timestamp']) + dyear = get_decimal_year(dt) + tm.append(dyear) + rate.append(float(row['value_60.0-120.0'])) + +plt.plot(tm, rate, color=cl, linewidth=1, linestyle='solid') + +#geminga_dt = datetime.date.fromisoformat('2023-04-16') +#day_of_year = geminga_dt.timetuple().tm_yday # returns 1 for January 1st +#geminga_tm=geminga_dt.year+(day_of_year/365) +#plt.axvline(x = geminga_tm, color = 'b', linewidth=3 , linestyle='dashed', label = 'Geminga scan') + +plt.ylabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal') +plt.xlabel('Year',fontsize=14, fontweight='normal') + +fout_png=input_file.replace("csv","png") +fig.savefig(products_dir+fout_png, bbox_inches='tight') +plt.show()