nustar/nuwavdet
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

base: nustar:master
Branches Tags
nustar:master
nustar:release_version
nustar:code_trimming
nustar:pre_github
..
compare: nustar:release_version
Branches Tags
nustar:master
nustar:release_version
nustar:code_trimming
nustar:pre_github

No commits in common. "master" and "release_version" have entirely different histories.
master ... release_ve

10 changed files with 653 additions and 826 deletions
Show Stats Expand all files Collapse all files

72
README.md
View File

@ -4,64 +4,34 @@ This pacakge is used to generate region masks separating any focused X-ray flux
## Installation ## Installation
This package is to be used with Python 3.x.x This package is to be used with Python 3.x.x
```bash
pip install git+http://heagit.cosmos.ru/nustar/nuwavdet.git
```
To update the package to the current version one should delete the previous version
```bash
pip uninstall nuwavdet
```
And simply repeat the intallation procedure again from the repository.
## Installation verification
If the installation was successful the package can be used with the following import:
```python ```python
from nuwavdet import nuwavdet as nw pip install git+https://github.com/andrey-rrousan/nuwavdet
```
To verify the installation we suggest running a simple script:
```python
from nuwavdet import nuwavdet as nw
print(nw.binary_array(2))
```
The output of the script should be
```bash
[[False False]
[False True]
[ True False]
[ True True]]
``` ```
## Main use ## Main use
The main functionality of the package is presented with a single function To use the package in your project, import it in by writing
```python ```python
nw.process(obs_path, thresh) from nuwavdet import nuwavdet as nw
``` ```
The main functionality of the pacakge is presented with a single function
```python
process(obs_path, thresh)
```
Inputs are string with path to the _cl.evt file to use and a tuple of thresholds, e.g. Inputs are string with path to the _cl.evt file to use and a tuple of thresholds, e.g.
```python ```python
nw.process('D:\\Data\\obs_cl.evt', (3, 2)) process('D:\\Data\\obs_cl.evt', (3, 2))
``` ```
The detailed script description of the data extraction with the script is presented in the examples folder of the repository. Outputs of the function are:
1. dictionary with some metadata and properties of the observation after mask generation procedure.
2. region array with mask in DET1 coordinate frame. Note that this mask is for numpy mask application so 1 corresponds to masked pixel and 0 otherwise.
3. custom bad pixel table with flagged pixels in RAW coordinates. It can be exported as fits file or each separate table can be acessed directly.
4. array with the sum of wavelet planes used in the processing.
The function nw.process returns severl python objects: Metadata about the observation file:
1. python-dictionary with some metadata and properties of the observation after mask generation procedure.
2. region array with mask in DET1 coordinate frame. Note that this mask is for numpy mask application so True (1) corresponds to masked pixel and False (0) otherwise.
3. custom bad pixel table with flagged pixels in RAW coordinates. It can be exported as fits file for further application to the nupipeline as fpma_userbpfile or fpmb_userbpfile.
4. array with the sum of wavelet planes for potential alternative applications.
Metadata about the observation returned by the nw.process is:
Observation metadata:
1. OBS_ID 1. OBS_ID
2. Detector 2. Detector
@ -77,8 +47,14 @@ Useful algorythm-related data:
## Other uses ## Other uses
Other possbile usecases are shown in the examples folder. You can process the cl.evt file by creating an Observation class object:
## Contact information ```python
obs = nw.Observation(path_to_evt_file)
```
If you have any questions or issues with the code, feel free to contact Andrey Mukhin: amukhin@cosmos.ru Additionally, the energy band in KeV to get events from can be passed as an argument. The default value is [3,20].
```python
obs = nuwavsource.Observation(path_to_evt_file,E_borders=[E_min,E_max])
```

54
examples/1_save_results.py
View File

@ -1,54 +0,0 @@
from nuwavdet import nuwavdet as nw
OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
THRESH = (3, 2)
SAVE_BADPIX_PATH = r'.//out//badpix.fits'
SAVE_REGION_PATH = r'.//out//region.fits'
SAVE_WAVSUM_PATH = r'.//out//wavsum.fits'
METADATA_PATH = r'.//out//metadata.csv'
METADATA_FITS_PATH = r'.//out//metadata.fits'
if __name__ == '__main__':
# PROCESS THE OBSERVATION WITH GIVEN THRESHOLD
result, region, region_raw, wav_sum = nw.process(OBS_PATH, thresh=THRESH)
# SAVE THE REGION BAD PIXEL FILES TO THE FITS FILE WITH NUPIPELINE
# COMPATIBLE FORMAT AND HEADERS.
region_raw.writeto(SAVE_BADPIX_PATH)
# SAVE REGION MASK AS A FITS IMAGE
nw.save_region(region, SAVE_REGION_PATH, overwrite=False)
# Note that the Python script uses numpy masked array with
# True (1) as as masked and False (0) as unmasked pixel.
# nw.save_region transfers the numpy masked array to
# conventional format with 1 for unmasked and 0 for masked pixel.
# However, if mask is used in the Python you need to transfer it back with
# numpy.logical_not(mask).
# SAVE WAVSUM ARRAY AS A FITS IMAGE
nw.fits.writeto(SAVE_WAVSUM_PATH, wav_sum, overwrite=False)
# SAVE METADATA
# WE SUGGEST SAVING ALL THE METADATA FOR SEVERAL OBSERVATIONS
# IN ONE FILE.
# CREATE CSV FILE TO SAVE DATA
# IF FILE ALREADY EXISTS YOU SHOULD REMOVE THIS BLOCK FROM YOUR CODE
table = {
'obs_id': [], 'detector': [], 'ra': [], 'dec': [],
'lon': [], 'lat': [], 't_start': [], 'exposure': [],
'count_rate': [], 'remaining_area': [], 'cash_stat': [],
'cash_stat_full': []
}
out_table = nw.DataFrame(table)
out_table.to_csv(METADATA_PATH)
# SAVE DATA TO CREATED CSV
nw.DataFrame(result, index=[0]).to_csv(METADATA_PATH, mode='a', header=False)
# TRANSFORM THE CSV TO FITS-TABLE
nw.csv_to_table(METADATA_PATH, METADATA_FITS_PATH)

33
examples/2_directory_processing.py
View File

@ -1,33 +0,0 @@
from nuwavdet import nuwavdet as nw
INPUT_FOLDER = r'path_to_directory'
OUTPUT_FOLDER = r'.//Output'
if __name__ == '__main__':
# BEGIN PROCESSING ALL THE OBSERVATIONS INSIDE THE FOLDER
nw.process_folder(input_folder=INPUT_FOLDER,
start_new_file='y',
fits_folder=OUTPUT_FOLDER,
thresh=(3, 2),
cpu_num=10
)
# IF THE PROCESSING WAS INTERRUPTED YOU CAN CONTINUE IT WITH THE SAME CODE
# BY CHANGING THE start_new_file TO 'n'.
# THE STRUCTURE OF THE FOLDER IS
# OUTPUT_FOLDER
# __overview.csv csv-table with observations metadata
# __overvies.fits fits-table with the same metadata
# __overview_skipped.csv csv-table with the skipped observations
# __Region folder for region mask images
# ____<obsid><DET>_region.fits
# __Region_raw folder for region masks in RAW coordinates
# ____<obsid><DET>_reg_raw.fits
# __Wav_sum folder for sum of wavelet layers
# ____<obsid><DET>_wav_sum.fits
# Note nw.process_folder uses multiprocessing with cpu_num cores.
# The number of cores can be manually chosen or automatically
# detected if cpu_num = 0.

23
examples/3_wavelet.py
View File

@ -1,23 +0,0 @@
from nuwavdet import nuwavdet as nw
OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
THRESH = (3, 2)
if __name__ == '__main__':
# CREATE THE OBSERVATION CLASS OBJECT
obs = nw.Observation(OBS_PATH)
# CALCULATE THE WAVLET LAYERS WITH GIVEN THRESHOLD
wav_layers = obs.wavdecomp(mode='atrous', occ_coeff=True, thresh=THRESH)
# ALL THE LAYERS CAN BE ACCESSED AS AN ELEMENT OF wav_layers VARIABLE
# wav_layers[0] for the 1st wavelet layer
# wav_layers[4] for 5th wavelet layer
# wav_layers[-1] for the last wavelet layer
# wav_layers[2:5] for the list of the layers from 3 to 5
# wav_layers[[1, 3, 5]] for the list of layers 2, 4 and 6
# To calculate the sum of wavelet layers one should use sum() method
# wav_layers[2:7].sum(0) returns a sum of layers from 3 to 7
# wav_layers[[1, 3, 5]].sum(0) returns a sum of layers 2, 4 and 6.

23
examples/4_cstat.py
View File

@ -1,23 +0,0 @@
from nuwavdet import nuwavdet as nw
import numpy as np
OBS_PATH = r'.//path_to_obs//nu<obsid><DET>01_cl.evt'
MASK_PATH = r'.//path_to_mask//<obsid><DET>.fits'
if __name__ == '__main__':
# CREATE THE OBSERVATION CLASS OBJECT
obs = nw.Observation(OBS_PATH)
# READ THE REGION MASK FILE
region = nw.fits.getdata(MASK_PATH)
# TRANSFORM REGION MASK DATA TO NUMPY MASK DATA (SEE 1_save_results.py).
region = np.logical_not(region.astype(bool))
# CREATE MASKED ARRAY CLASS OBJECT
masked_data = np.ma.masked_array(obs, mask=region)
# CALCULATE THE CSTAT ON THE MASKED DATA
print(nw.сstat(masked_data.mean(), masked_data))

BIN
nuwavdet/__pycache__/__init__.cpython-39.pyc Normal file
View File

Binary file not shown.

BIN
nuwavdet/__pycache__/nuwavdet.cpython-39.pyc Normal file
View File

Binary file not shown.

BIN
nuwavdet/__pycache__/nuwavsource.cpython-39.pyc Normal file
View File

Binary file not shown.

66
nuwavdet/nuwavdet.py
View File

@ -30,14 +30,6 @@ def get_link_list(folder: str, sort_list: bool = True) -> list[str]:
return np.array(links) return np.array(links)
def csv_to_table(csv_path, fits_path):
"""
Transform the csv table to fits table with astropy.
"""
csv_file = read_csv(csv_path, index_col=0, dtype={'obs_id': str})
Table.from_pandas(csv_file).write(fits_path, overwrite=True)
def binary_array(num: int) -> list[list[bool]]: def binary_array(num: int) -> list[list[bool]]:
""" """
Returns list of all possible combinations of num of bool values. Returns list of all possible combinations of num of bool values.
@ -207,8 +199,7 @@ def count_binning(array, count_per_bin: int = 2):
def cstat(expected, data: list, count_per_bin: int = 2) -> float: def cstat(expected, data: list, count_per_bin: int = 2) -> float:
_data = data.flatten() _data = data.flatten()
if type(data) is np.ma.masked_array: _data = _data[_data.mask == False]
_data = _data[_data.mask == False]
_expected = expected _expected = expected
c_stat = 0 c_stat = 0
bin_sum_array, bin_count_array = count_binning(_data, count_per_bin) bin_sum_array, bin_count_array = count_binning(_data, count_per_bin)
@ -244,7 +235,7 @@ class Observation:
resized_coeff = (coeff).reshape(2, 2).repeat(180, 0).repeat(180, 1) resized_coeff = (coeff).reshape(2, 2).repeat(180, 0).repeat(180, 1)
return resized_coeff return resized_coeff
def get_data(self, file, E_borders=[3, 20], generate_mask=True): def get_data(self, file, E_borders=[3, 20]):
""" """
Returns masked array with DET1 image data for given energy band. Returns masked array with DET1 image data for given energy band.
Mask is created from observations badpix tables and to mask the border and gaps. Mask is created from observations badpix tables and to mask the border and gaps.
@ -257,12 +248,10 @@ class Observation:
data_mask = data[np.logical_not(idx_mask)] data_mask = data[np.logical_not(idx_mask)]
build_hist = lambda array: np.histogram2d(array['DET1Y'], array['DET1X'], 360, [[0, 360], [0, 360]])[0] build_hist = lambda array: np.histogram2d(array['DET1Y'], array['DET1X'], 360, [[0, 360], [0, 360]])[0]
output = build_hist(data_output) output = build_hist(data_output)
if generate_mask: mask = build_hist(data_mask)
mask = build_hist(data_mask) mask = np.logical_or(mask, add_borders(output))
mask = np.logical_or(mask, add_borders(output)) mask = np.logical_or(mask, self.get_bad_pix(file))
mask = np.logical_or(mask, self.get_bad_pix(file)) return output, mask
return output, mask
return output
def get_bad_pix(self, file, threshold=0.9): def get_bad_pix(self, file, threshold=0.9):
""" """
@ -287,7 +276,7 @@ class Observation:
correction_poiss = np.random.poisson(corr*array, corr.shape) correction_poiss = np.random.poisson(corr*array, corr.shape)
return array + correction_poiss return array + correction_poiss
def wavdecomp(self, mode='gauss', thresh=0, occ_coeff=False): def wavdecomp(self, mode='gauss', thresh=False, occ_coeff=False):
""" """
Performs a wavelet decomposition of image. Performs a wavelet decomposition of image.
""" """
@ -341,7 +330,7 @@ class Observation:
""" """
Returns a hdu_list with positions of masked pixels in RAW coordinates. Returns a hdu_list with positions of masked pixels in RAW coordinates.
""" """
y_region, x_region = np.where(region) x_region, y_region = np.where(region)
hdus = [] hdus = []
for i in range(4): for i in range(4):
current_dir = os.path.dirname(__file__) current_dir = os.path.dirname(__file__)
@ -401,6 +390,7 @@ def process(obs_path, thresh):
Arguments: path to the file of interest and threshold, Arguments: path to the file of interest and threshold,
e.g. process('D:\\Data\\obs_cl.evt', (3, 2)) e.g. process('D:\\Data\\obs_cl.evt', (3, 2))
""" """
bin_num = 6
table = { table = {
'obs_id': [], 'detector': [], 'ra': [], 'dec': [], 'obs_id': [], 'detector': [], 'ra': [], 'dec': [],
@ -415,19 +405,18 @@ def process(obs_path, thresh):
dec=obs.dec*u.deg, dec=obs.dec*u.deg,
frame='fk5').transform_to('galactic') frame='fk5').transform_to('galactic')
lon, lat = sky_coord.l.value, sky_coord.b.value lon, lat = sky_coord.l.value, sky_coord.b.value
useful_bin_num = 6 rem_signal, rem_area, poiss_comp, rms = np.zeros((4, 2**bin_num))
rem_signal, rem_area, poiss_comp, rms = np.zeros((4, 2**useful_bin_num))
region = np.zeros(obs.data.shape, dtype=bool) region = np.zeros(obs.data.shape, dtype=bool)
region_raw = -1 region_raw = -1
rem_region = np.logical_and(region, np.logical_not(obs.data.mask)) rem_region = np.logical_and(region, np.logical_not(obs.data.mask))
masked_obs = np.ma.masked_array(obs.data, mask=region) masked_obs = np.ma.masked_array(obs.data, mask=region)
good_lvl = np.zeros(useful_bin_num, dtype=bool) good_lvl = np.zeros(bin_num, dtype=bool)
good_idx = 0
if obs.exposure > 1000: if obs.exposure > 1000:
wav_obs = obs.wavdecomp('atrous', thresh, occ_coeff=True) wav_obs = obs.wavdecomp('atrous', thresh, occ_coeff=True)
wav_sum = wav_obs[2:-1].sum(0) wav_sum = wav_obs[2:-1].sum(0)
occ_coeff = obs.get_coeff() occ_coeff = obs.get_coeff()
binary_arr = binary_array(useful_bin_num) binary_arr = binary_array(bin_num)
good_idx = len(binary_arr) - 1
for idx, lvl in enumerate(binary_arr): for idx, lvl in enumerate(binary_arr):
try: try:
@ -445,9 +434,12 @@ def process(obs_path, thresh):
rms[idx] = np.sqrt(((masked_obs-masked_obs.mean())**2).mean()) rms[idx] = np.sqrt(((masked_obs-masked_obs.mean())**2).mean())
for idx in range(len(poiss_comp)): for idx in range(len(poiss_comp)):
if ((poiss_comp[idx] < poiss_comp[-1] + 0.05) and if ((poiss_comp[idx] < poiss_comp[good_idx]) and
(rem_area[idx] > rem_area[good_idx])): (poiss_comp[idx] < poiss_comp[-1] + 0.05) and
(rem_area[idx] > rem_area[-1])):
good_idx = idx good_idx = idx
if good_idx == 0:
good_idx = len(binary_arr) - 1
good_lvl = binary_arr[good_idx] good_lvl = binary_arr[good_idx]
try: try:
@ -491,7 +483,7 @@ def process(obs_path, thresh):
] ]
for key, value in zip(table.keys(), to_table): for key, value in zip(table.keys(), to_table):
table[key] = value table[key] = [value]
return table, region.astype(int), region_raw, wav_sum return table, region.astype(int), region_raw, wav_sum
except TypeError: except TypeError:
return obs_path, -1, -1, -1 return obs_path, -1, -1, -1
@ -502,7 +494,7 @@ def _process_multi(args):
def process_folder(input_folder=None, start_new_file=None, fits_folder=None, def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
thresh=None, cpu_num=0): thresh=None):
""" """
Generates a fits-table of parameters, folder with mask images in DET1 and Generates a fits-table of parameters, folder with mask images in DET1 and
BADPIX tables in RAW for all observations in given folder. BADPIX tables in RAW for all observations in given folder.
@ -594,14 +586,6 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
# START PROCESSING # START PROCESSING
print('Started processing...') print('Started processing...')
num = 0 num = 0
if cpu_num == 0:
cpu_num = cpu_count()
elif cpu_num < 0:
raise ValueError('cpu_num must be a positive integer')
elif cpu_num > cpu_count():
print('Chosen cpu_num exceed the number of CPU cores. Using cpu_count() instead.')
cpu_num = cpu_count()
for group_idx in range(len(obs_list)//group_size+1): for group_idx in range(len(obs_list)//group_size+1):
print(f'Started group {group_idx}') print(f'Started group {group_idx}')
group_list = obs_list[group_size*group_idx:min(group_size*(group_idx+1), len(obs_list))] group_list = obs_list[group_size*group_idx:min(group_size*(group_idx+1), len(obs_list))]
@ -609,7 +593,7 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
os.stat(file).st_size/2**20 os.stat(file).st_size/2**20
for file in group_list for file in group_list
]).max() ]).max()
process_num = (cpu_num if max_size < 50 else (cpu_num//2 if max_size < 200 else cpu_num//4 if max_size < 1000 else cpu_num//8)) process_num = (cpu_count() if max_size < 50 else (cpu_count()//2 if max_size < 200 else (cpu_count()//4 if max_size < 1000 else 1)))
print(f"Max file size in group is {max_size:.2f}Mb, create {process_num} processes") print(f"Max file size in group is {max_size:.2f}Mb, create {process_num} processes")
with get_context('spawn').Pool(processes=process_num) as pool: with get_context('spawn').Pool(processes=process_num) as pool:
packed_args = map(lambda _: (_, thresh), group_list) packed_args = map(lambda _: (_, thresh), group_list)
@ -620,14 +604,14 @@ def process_folder(input_folder=None, start_new_file=None, fits_folder=None,
num += 1 num += 1
continue continue
obs_name = str(result['obs_id'])+result['detector'] obs_name = str(result['obs_id'][0])+result['detector'][0]
if result['exposure'] < 1000: if result['exposure'][0] < 1000:
print(f'{num:>3} {obs_name} is skipped. Exposure < 1000') print(f'{num:>3} {obs_name} is skipped. Exposure < 1000')
DataFrame(result, index=[0]).to_csv(os.path.join(fits_folder, 'overview_skipped.csv'), mode='a', header=False) DataFrame(result).to_csv(os.path.join(fits_folder, 'overview_skipped.csv'), mode='a', header=False)
num += 1 num += 1
continue continue
DataFrame(result, index=[0]).to_csv(os.path.join(fits_folder, 'overview.csv'), mode='a', header=False) DataFrame(result).to_csv(os.path.join(fits_folder, 'overview.csv'), mode='a', header=False)
save_region(region, os.path.join(region_folder, f'{obs_name}_region.fits'), overwrite=True) save_region(region, os.path.join(region_folder, f'{obs_name}_region.fits'), overwrite=True)
region_raw.writeto(os.path.join(region_raw_folder, f'{obs_name}_reg_raw.fits'), overwrite=True) region_raw.writeto(os.path.join(region_raw_folder, f'{obs_name}_reg_raw.fits'), overwrite=True)
fits.writeto(os.path.join(wav_sum_folder, f'{obs_name}_wav_sum.fits'), wav_sum, overwrite=True) fits.writeto(os.path.join(wav_sum_folder, f'{obs_name}_wav_sum.fits'), wav_sum, overwrite=True)

6
setup.py
View File

@ -5,10 +5,10 @@ with open("README.md", "r") as fh:
setuptools.setup( setuptools.setup(
name="nuwavdet", name="nuwavdet",
version="0.1.1", version="0.1.0",
author="Andrey Mukhin", author="Andrey Mukhin",
author_email="amukhin@cosmos.ru", author_email="amukhin@phystech.edu",
description="A package for source exclusion in NuSTAR observation data using wavelet decomposition", description="A package for source exclusion in NuStar observation data using wavelet decomposition",
long_description=long_description, long_description=long_description,
long_description_content_type="text/markdown", long_description_content_type="text/markdown",
url="https://github.com/andrey-rrousan/nuwavdet", url="https://github.com/andrey-rrousan/nuwavdet",