Commit de57eedd authored by Carl Schaffer's avatar Carl Schaffer
Browse files

Merge branch 'reformatting' into 'master'

Reformatting and Cleaning

See merge request !140
parents f5c0a833 17de514c
......@@ -4,4 +4,4 @@
"""
Created by schaffer at 7/30/19
"""
\ No newline at end of file
"""
......@@ -42,11 +42,11 @@ class OrderedLoader(yaml.Loader):
def get_separators(yaml_file):
try:
f = open(yaml_file, 'rb')
f = open(yaml_file, "rb")
yamldata = f.read()
f.close()
except:
print('Could not open file {:}'.format(yaml_file))
print("Could not open file {:}".format(yaml_file))
return None
sep_data = yaml.load(yamldata, OrderedLoader)["separators"]
......@@ -62,31 +62,31 @@ def get_generic_header(filename=settings.header_template, version="reconstructio
"""
try:
f = open(filename, 'rb')
f = open(filename, "rb")
yamldata = f.read()
f.close()
except:
print('Could not open file {:}'.format(filename))
print("Could not open file {:}".format(filename))
return None
headerdata = yaml.load(yamldata, OrderedLoader)
try:
header = fits.Header(headerdata[version].keys())
except:
print('Version {:} not found'.format(version))
print("Version {:} not found".format(version))
return None
headerversion = headerdata[version]
for card in headerversion:
try:
gencard = headerversion.get(card)
standard = gencard['DEFAULT']
comment = gencard.get('COMMENT', '')
standard = gencard["DEFAULT"]
comment = gencard.get("COMMENT", "")
if comment is None:
comment = ''
comment = ""
header[card] = (standard, comment)
except:
print('Could not update header')
print("Could not update header")
print(card, gencard)
return header
......@@ -99,11 +99,11 @@ def get_generic_patches(filename, version):
:return: patches ordered dict
"""
try:
f = open(filename, 'rb')
f = open(filename, "rb")
yamldata = f.read()
f.close()
except:
print('Could not open file {:}'.format(filename))
print("Could not open file {:}".format(filename))
return None
patchesdata = yaml.load(yamldata, OrderedLoader)
......@@ -111,7 +111,7 @@ def get_generic_patches(filename, version):
try:
patches = patchesdata[version]
except:
print('Version {:} not found'.format(version))
print("Version {:} not found".format(version))
return None
return patches
......@@ -133,7 +133,7 @@ def apply_patches(new_header, old_header, patches):
try:
new_header = apply_patch(new_header, old_header, patch)
except:
print('Could not apply patch {0}'.format(patch))
print("Could not apply patch {0}".format(patch))
pass
return new_header
......@@ -158,21 +158,21 @@ def apply_patch(new_header, old_header, patch):
"""
try:
kind = patch[1]['KIND']
old = patch[1]['OLD']
new = patch[1]['NEW']
function = patch[1]['FUNCTION']
if kind == 'UPDATEVALUE':
function = function.replace('OLD', "'" + old + "'")
function = function.replace('NEW', "'" + new + "'")
kind = patch[1]["KIND"]
old = patch[1]["OLD"]
new = patch[1]["NEW"]
function = patch[1]["FUNCTION"]
if kind == "UPDATEVALUE":
function = function.replace("OLD", "'" + old + "'")
function = function.replace("NEW", "'" + new + "'")
else:
print('Unknown Kind: {0}'.format(kind))
print("Unknown Kind: {0}".format(kind))
exec(function)
except KeyError:
print(f"Error Executing {function}")
except:
print('Could not evaluate patch: {0}'.format(function))
print("Could not evaluate patch: {0}".format(function))
pass
return new_header
......@@ -219,7 +219,7 @@ def insert_wcs_kws(header, filename):
# add wcs kws
x, y = get_coords(filename)
data = fits.getdata(filename)
center = [int(ax / 2.) for ax in data.shape]
center = [int(ax / 2.0) for ax in data.shape]
center.reverse() # fix x and y order
npix_x, npix_y = data.shape[1], data.shape[0] # fix x and y order
......@@ -234,13 +234,17 @@ def insert_wcs_kws(header, filename):
CRPIX1=center[0],
CRVAL1=float(x),
CDELT1=step_width_x,
CSYER1=fov_x * .25 * step_width_x, # assume coordinate error of 25 percent of image
CSYER1=fov_x
* 0.25
* step_width_x, # assume coordinate error of 25 percent of image
CTYPE2="HPLT-TAN",
CUNIT2="arcsec",
CRPIX2=center[1],
CRVAL2=float(y),
CDELT2=step_width_y,
CSYER2=fov_y * .25 * step_width_y, # assume coordinate error of 25 percent of image
CSYER2=fov_y
* 0.25
* step_width_y, # assume coordinate error of 25 percent of image
# Compatibility for SunPy Maps
HGLN_OBS=0,
HGLT_OBS=0,
......@@ -268,17 +272,17 @@ def insert_separators(new_header):
for sep in sep_list:
for sep_contents, sep_config in sep.items():
target = sep_config["after"]
new_header.add_comment('-' * 70, after=target)
new_header.add_comment("-" * 70, after=target)
if "additional_content" in sep_config:
content = sep_config['additional_content']
content = sep_config["additional_content"]
lines = content.split('\n')
lines = content.split("\n")
lines.reverse()
for line in lines:
new_header.add_comment(line, after=target)
new_header.add_comment(sep_contents.center(70, '-'), after=target)
new_header.add_comment('-' * 70, after=target)
new_header.add_comment(sep_contents.center(70, "-"), after=target)
new_header.add_comment("-" * 70, after=target)
def convert_bbi_header(filename):
......@@ -291,15 +295,15 @@ def convert_bbi_header(filename):
"""
header = fits.getheader(filename)
bbi_header = get_generic_header(settings.header_template, 'reconstruction')
bbi_header = get_generic_header(settings.header_template, "reconstruction")
bbi_fn_template = settings.fn_template
patches = get_generic_patches(settings.header_patches_template, 'reconstruction')
patches = get_generic_patches(settings.header_patches_template, "reconstruction")
# Check how much changes are necessary and inform user
diff = fits.HeaderDiff(header, bbi_header)
if diff.identical:
print('Header is identical to generic header')
print("Header is identical to generic header")
return header
new_header = bbi_header.copy()
......@@ -315,21 +319,23 @@ def convert_bbi_header(filename):
new_header = apply_patches(new_header, header, patches)
date_section = re.search(r"\d{8}[_\-]\d{6}", basename(filename)).group(0)
timestamp = datetime.datetime.strptime(date_section, '%Y%m%d-%H%M%S')
wavelength = new_header['WAVELNTH']
timestamp = datetime.datetime.strptime(date_section, "%Y%m%d-%H%M%S")
wavelength = new_header["WAVELNTH"]
pointing_id = get_pointing_id(filename)
new_header["POINT_ID"] = pointing_id
new_header["EXTNAME"] = "_".join([new_header["POINT_ID"], f"WL{wavelength:03.0f}"])
new_header["FILTER"] = get_filter(wavelength)
file_mode = "reco_bb" if (('Zyla-1' in filename) or ('narrow' in filename)) else "reco_nb"
file_mode = (
"reco_bb" if (("Zyla-1" in filename) or ("narrow" in filename)) else "reco_nb"
)
fits_filename = bbi_fn_template.format(
timestamp.strftime("%Y%m%d-%H%M%S"),
pointing_id.split("_")[-1],
int(wavelength),
file_mode,
"fits"
"fits",
)
new_filename = fits_filename
......@@ -341,13 +347,15 @@ def convert_bbi_header(filename):
raise KeyError
# speck = pipeline.speckle(filename + '.log')
except:
print('Could not find metadata')
print("Could not find metadata")
speck = None
if speck is not None:
estimated_alpha = speck.get_sub_alpha()
estimated_alpha = estimated_alpha.mean((0, 1))[0]
estimated_fried = estimated_alpha * float(speck.parameter['telescopediameter_mm']) / 10.
new_header['ATMOS_R0'] = float(f"{estimated_fried:03.3f}")
estimated_fried = (
estimated_alpha * float(speck.parameter["telescopediameter_mm"]) / 10.0
)
new_header["ATMOS_R0"] = float(f"{estimated_fried:03.3f}")
insert_wcs_kws(new_header, filename)
add_history(new_header, __file__, "header_fix")
......@@ -355,7 +363,7 @@ def convert_bbi_header(filename):
return new_header
if __name__ == '__main__':
if __name__ == "__main__":
infile = "/archive/bbi/new/20170519/Zyla1/RC/RC_data_andorZyla-1_20170519-123314144_HAlpha_656.fits"
header_template = "/dat/schaffer/projects/bbi_header/generic_bbi_header.yaml"
outpath = "/tmp/"
......@@ -370,7 +378,7 @@ def run_bbi_header_conversion(args):
"--folder_out",
dest="folder_out",
default="/dat/sdc/bbi/",
help="Output folder, all output files will be stored here, (default: %(default)s)"
help="Output folder, all output files will be stored here, (default: %(default)s)",
)
args = parser.parse_args(args)
in_folder = args.folder_in
......@@ -408,9 +416,16 @@ def run_bbi_header_conversion(args):
out_dir = join(out_folder, new_header["POINT_ID"] + f"_WL{wl}")
setup_outdir(out_dir)
out_file = join(out_dir, new_header['FILENAME'])
out_file = join(out_dir, new_header["FILENAME"])
data = fits.getdata(infile)
flipped = np.flip(data, axis=0)
fits.writeto(out_file, flipped, new_header, overwrite=True, checksum=True, output_verify='warn')
fits.writeto(
out_file,
flipped,
new_header,
overwrite=True,
checksum=True,
output_verify="warn",
)
......@@ -13,7 +13,7 @@ import tqdm
from kis_tools.generic.settings.settings import get_settings
from kis_tools.util.util import date_from_fn
settings = get_settings('bbi')
settings = get_settings("bbi")
path_obs_csv = settings.path_obs_csv
......@@ -31,7 +31,7 @@ def copy_file(infile, run_number):
def get_obs_info(filename):
date = date_from_fn(filename)
with open(path_obs_csv, 'r') as csv_file:
with open(path_obs_csv, "r") as csv_file:
obs_csv = csv.reader(csv_file)
header = next(obs_csv)
pos_date_end = header.index("DATE-END")
......@@ -59,7 +59,7 @@ def get_obs_info(filename):
def make_run_folders():
path = join(dirname(__file__), "bbi_observations_with_coordinates.csv"),
path = (join(dirname(__file__), "bbi_observations_with_coordinates.csv"),)
df = pd.read_csv(
path,
......
......@@ -24,9 +24,7 @@ def get_times(folder):
""" extract dates from all files from a given folder"""
files = glob(join(folder, OBS_FILE_PATTERN))
stamps = sorted([f.split("_")[2] for f in files])
times = np.array(
[datetime.datetime.strptime(s, "%Y%m%d-%H%M%S%f") for s in stamps]
)
times = np.array([datetime.datetime.strptime(s, "%Y%m%d-%H%M%S%f") for s in stamps])
return times
......@@ -78,9 +76,7 @@ if __name__ == "__main__":
beg, end = dates[0], dates[1]
name = beg.strftime("%Y%m%d")
name += f"_{i:03d}"
observations.append(
{"POINT_ID": name, "DATE-END": end, "DATE-BEG": beg}
)
observations.append({"POINT_ID": name, "DATE-END": end, "DATE-BEG": beg})
df = pd.DataFrame(observations)
df = df[["POINT_ID", "DATE-BEG", "DATE-END"]]
df["DURATION"] = df["DATE-END"] - df["DATE-BEG"]
......
......@@ -20,7 +20,7 @@ from tqdm import tqdm
from .util import get_wl
from ..util.util import merge_to_gif, gif_to_mp4
mpl.use('Agg')
mpl.use("Agg")
import matplotlib.pyplot as plt # noqa: E402
from ..generic.fits import FitsFile # noqa: E402
from sunpy.map import Map # noqa: E402
......@@ -31,21 +31,23 @@ import numpy as np # noqa: E402
def loc_preview(fits_file, fn_out):
"""Generate a preview showing the center of the observation on a up to date image of the soalr disk."""
m = Map(fits_file)
x, y = (m.top_right_coord.Tx + m.bottom_left_coord.Tx) / 2, (m.top_right_coord.Ty + m.bottom_left_coord.Ty) / 2
x, y = (m.top_right_coord.Tx + m.bottom_left_coord.Tx) / 2, (
m.top_right_coord.Ty + m.bottom_left_coord.Ty
) / 2
time = FitsFile(fits_file).obs_time
with catch_warnings():
warnings.simplefilter("ignore")
fig, ax = make_loc_plot(x.value, y.value, time)
# plt.tight_layout()
fig.savefig(fn_out, bbox_inches="tight", pad_inches=-.01, dpi=200)
fig.savefig(fn_out, bbox_inches="tight", pad_inches=-0.01, dpi=200)
return fn_out
def map_preview(fits_file, fn_out):
"""Generate a preview of the reconstructed data."""
m = Map(fits_file)
m.plot_settings['vmin'] = np.percentile(m.data, 0.1)
m.plot_settings['vmax'] = np.percentile(m.data, 99.9)
m.plot_settings["vmin"] = np.percentile(m.data, 0.1)
m.plot_settings["vmax"] = np.percentile(m.data, 99.9)
_ = m.plot()
plt.gcf().savefig(fn_out, bbox_inches="tight", pad_inches=0.1, dpi=200)
......@@ -56,10 +58,11 @@ def alpha_preview(fits_file, fn_out):
"""Generate a preview of the estimated seeing quality over the different reconstruction subfields."""
datetime = FitsFile(fits_file).obs_time
format_keys = dict(
date=datetime.strftime("%Y%m%d"),
time=datetime.strftime("%H%M%S")
date=datetime.strftime("%Y%m%d"), time=datetime.strftime("%H%M%S")
)
log = glob(
"/archive/bbi/new/{date}/Zyla1/RC/RC*{date}-{time}*.log".format(**format_keys)
)
log = glob("/archive/bbi/new/{date}/Zyla1/RC/RC*{date}-{time}*.log".format(**format_keys))
if not log:
raise IOError(f"Could not find reconstruction results for {fits_file}")
speck = pipeline.speckle(log[0])
......@@ -70,8 +73,7 @@ def alpha_preview(fits_file, fn_out):
def bbi_gen_gif(folder):
"""Check a folder and merge all images of the same type into gifs and mp4s"""
for pattern in ['reco', 'location', 'alpha']:
for pattern in ["reco", "location", "alpha"]:
files = glob(join(folder, f"*{pattern}*.png"))
if not files:
continue
......@@ -82,16 +84,17 @@ def bbi_gen_gif(folder):
for group, grouped_files in grouped:
grouped_files = list(grouped_files)
fn = basename(grouped_files[0])
fn = re.sub(r'-\d{6}', '', fn)
fn = fn.replace('.png', '.gif')
fn = re.sub(r"-\d{6}", "", fn)
fn = fn.replace(".png", ".gif")
print(fn)
merge_to_gif(sorted(grouped_files), join(folder, fn))
gif_to_mp4(join(folder, fn))
def run_bbi_preview(cmd_line_args):
mpl.use('Agg')
mpl.use("Agg")
from ..interface.commandline_interface import add_arg_infiles
parser = argparse.ArgumentParser()
add_arg_infiles(parser)
args = parser.parse_args(cmd_line_args)
......@@ -107,7 +110,7 @@ def run_bbi_preview(cmd_line_args):
map_preview(infile, fn_out_reco)
loc_preview(infile, fn_out_loc)
alpha_preview(infile, fn_out_alpha)
plt.close('all')
plt.close("all")
# folder = dirname(infiles[0])
# bbi_gen_gif(folder)
\ No newline at end of file
# bbi_gen_gif(folder)
......@@ -3,4 +3,4 @@ import re
def get_wl(fn):
"""Extract wavelength from string"""
return int(re.search(r'WL(\d+)', fn).group(1))
\ No newline at end of file
return int(re.search(r"WL(\d+)", fn).group(1))
......@@ -58,9 +58,9 @@ class ChrotelFitsFile(FitsFile):
self.filename = abspath(filename)
self.basename = basename(filename)
self.dir = dirname(self.filename)
self.preview_dir = ''
self.preview_dir = ""
self.has_preview_dir = False
self.movie_dir = ''
self.movie_dir = ""
self.has_movie_dir = False
self.parse_filepath()
......@@ -69,7 +69,7 @@ class ChrotelFitsFile(FitsFile):
def parse_filepath(self):
sub = self.basename.split('_')
sub = self.basename.split("_")
self.date = date_from_fn(self.basename)
self.filt = sub[1]
......@@ -81,38 +81,55 @@ class ChrotelFitsFile(FitsFile):
def get_dir_struct(self):
"""
Retrieve directory structure.
"""
self.preview_dir = abspath(join(self.dir, "../../../../", "jpeg", self.dir_date))
self.movie_dir = abspath(join(self.dir, "../../../../../", "movie", self.dir_date))
self.preview_dir = abspath(
join(self.dir, "../../../../", "jpeg", self.dir_date)
)
self.movie_dir = abspath(
join(self.dir, "../../../../../", "movie", self.dir_date)
)
self.has_preview_dir = isdir(self.preview_dir)
self.has_movie_dir = isdir(self.movie_dir)
@property
def preview_file(self):
file_path = ''
file_path = ""
if self.has_preview_dir:
file_path = join(self.preview_dir, self.basename.replace("fits.gz", "jpg"))
if not isfile(file_path):
file_path = ''
file_path = ""
return file_path
@property
def movie_file(self):
file_path = ''
file_path = ""
if self.has_movie_dir:
sub = self.basename.split('_')
file_name = '_'.join([sub[0], sub[1], sub[3]]) + '.avi'
sub = self.basename.split("_")
file_name = "_".join([sub[0], sub[1], sub[3]]) + ".avi"
file_path = join(self.movie_dir, file_name)
if not isfile(file_path):
file_path = ''
file_path = ""
return file_path
@property
def center(self):
def get_axis(axis_number):
ref_pixels = [int(self.header[k]) for k in self.header.keys() if f"CRPIX{axis_number}" in k]
ref_values = [float(self.header[k]) for k in self.header.keys() if f"CRVAL{axis_number}" in k]
step_values = [float(self.header[k]) for k in self.header.keys() if f"CDELT{axis_number}" in k]
ref_pixels = [
int(self.header[k])
for k in self.header.keys()
if f"CRPIX{axis_number}" in k
]
ref_values = [
float(self.header[k])
for k in self.header.keys()
if f"CRVAL{axis_number}" in k
]
step_values = [
float(self.header[k])
for k in self.header.keys()
if f"CDELT{axis_number}" in k
]
extent = int(self.query(f"NAXIS{axis_number}"))
assert len(ref_pixels) == len(ref_values)
......@@ -135,11 +152,22 @@ class ChrotelFitsFile(FitsFile):
@property
def box(self):
def get_minmax_axis(axis_number):
ref_pixels = [int(self.header[k]) for k in self.header.keys() if f"CRPIX{axis_number}" in k]
ref_values = [float(self.header[k]) for k in self.header.keys() if f"CRVAL{axis_number}" in k]
step_values = [float(self.header[k]) for k in self.header.keys() if f"CDELT{axis_number}" in k]
ref_pixels = [
int(self.header[k])
for k in self.header.keys()
if f"CRPIX{axis_number}" in k
]
ref_values = [
float(self.header[k])
for k in self.header.keys()
if f"CRVAL{axis_number}" in k
]
step_values = [
float(self.header[k])
for k in self.header.keys()
if f"CDELT{axis_number}" in k
]
extent = int(self.query(f"NAXIS{axis_number}"))
assert len(ref_pixels) == len(ref_values)
......
from .fits import FitsFile
\ No newline at end of file
from .fits import FitsFile
......@@ -42,7 +42,6 @@ class FitsFile:
except KeyError as e:
raise AttributeError
def _get_value(self, keywords):
"""
Try to retrieve a value from the header, allows multiple keys to be passed as a list. All given
......@@ -137,7 +136,9 @@ class FitsFile:
def _calc_heliocentric_angles(self):
coords = self.sunpy_map.center
mu, theta = calc_heliocentric_angle(coords.Tx.arcsec, coords.Ty.arcsec, date=self.obs_time)
mu, theta = calc_heliocentric_angle(
coords.Tx.arcsec, coords.Ty.arcsec, date=self.obs_time
)
self._mu = mu
self._theta = theta
......@@ -161,7 +162,7 @@ class FitsFile:
obs_time: observation time given in filename
"""
time_string = re.search(r"(\d{8}[\-_]\d{6})", self.filename).group(1)
time_string = time_string.replace('-', '_')
time_string = time_string.replace("-", "_")
obs_time = datetime.datetime.strptime(time_string, "%Y%m%d_%H%M%S")
return obs_time
......
......@@ -58,16 +58,19 @@ class ChroTelSettings(SDCSettings, metaclass=MetaClassSettings):
class BBISettings(SDCSettings, metaclass=MetaClassSettings):
def __init__(self):
self.settings_file = join(_BASE_DIR, "generic/settings/db_structure.yaml")
super().__init__(self.settings_file)
# extract observations config from yaml
self.path_obs_csv = join(_BASE_DIR,'bbi', self.obs_csv)
self.path_obs_csv = join(_BASE_DIR, "bbi", self.obs_csv)
self.header_template = join(self._source_directory, self.header_template_fn)
self.header_patches_template = join(self._source_directory, self.header_patches_fn)
self.observations = pd.read_csv(self.path_obs_csv, parse_dates=["DATE-BEG", "DATE-END"])
self.header_patches_template = join(
self._source_directory, self.header_patches_fn
)
self.observations = pd.read_csv(
self.path_obs_csv, parse_dates=["DATE-BEG", "DATE-END"]
)
def get_settings(instrument_name: str):
......@@ -77,7 +80,6 @@ def get_settings(instrument_name: str):
"lars": SDCSettings,
"gris": SDCSettings,
"chrotel": ChroTelSettings,
}
if instrument_name in settings_class:
......@@ -86,4 +88,6 @@ def get_settings(instrument_name: str):
settings_instance = settings_class()
return settings_instance