artxc/monitor
1
0
Fork 0
generated from erosita/uds
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

final

Browse Source
This commit is contained in:
Roman Krivonos
2025-10-09 17:37:20 +03:00
parent 879d09ab65
commit 7b9ab710e7
6 changed files with 562 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
README.md
View File

@@ -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 ```
## Работа с данными

3
monitor/monitor/config.py
View File

@@ -0,0 +1,3 @@
threshold=1.1
download_link='https://monitor.srg.cosmos.ru/download?type=orig&range=60.0-120.0'

161
monitor/monitor/utils.py
View File

@@ -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')

77
scripts/download.py Executable file
View File

@@ -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}")

205
scripts/gti.py Executable file
View File

@@ -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)

111
scripts/plot.py Executable file
View File

@@ -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()