192 lines
7.3 KiB
Python
192 lines
7.3 KiB
Python
import os
|
|
from collections import defaultdict
|
|
from multiprocessing import Queue, Process
|
|
import re
|
|
import traceback
|
|
|
|
import tqdm
|
|
import numpy as np
|
|
from matplotlib import pyplot as plt
|
|
from astropy.io import fits
|
|
from skimage.transform import hough_circle, hough_circle_peaks
|
|
from skimage.feature import canny
|
|
from skimage.filters import gaussian
|
|
from skimage.morphology import skeletonize
|
|
import cv2 as cv
|
|
|
|
def lowpriority():
|
|
""" Set the priority of the process to below-normal."""
|
|
|
|
import sys
|
|
try:
|
|
sys.getwindowsversion()
|
|
except AttributeError:
|
|
isWindows = False
|
|
else:
|
|
isWindows = True
|
|
|
|
if isWindows:
|
|
# Based on:
|
|
# "Recipe 496767: Set Process Priority In Windows" on ActiveState
|
|
# http://code.activestate.com/recipes/496767/
|
|
import win32api,win32process,win32con # pywin32
|
|
|
|
pid = win32api.GetCurrentProcessId()
|
|
handle = win32api.OpenProcess(win32con.PROCESS_ALL_ACCESS, True, pid)
|
|
win32process.SetPriorityClass(handle, win32process.BELOW_NORMAL_PRIORITY_CLASS)
|
|
else:
|
|
import os
|
|
|
|
os.nice(1)
|
|
|
|
measurement_names = ["094", "131", "171", "195", "284", "304"]
|
|
thresholds = [0.050, 0.10 , 1.00, 1.40, 0.90, 2.50]
|
|
max_center_skew = 7
|
|
circle_hough_thresh = 0.7
|
|
expected_dims = 1280
|
|
ratio_above_thresh_max = 0.5
|
|
|
|
def filter_fits(work_queue):
|
|
lowpriority()
|
|
while True:
|
|
try:
|
|
job = work_queue.get()
|
|
if job is None:
|
|
return
|
|
measurement = job.split("dr_suvi-l2-")[1].split("_")[0][2:]
|
|
idx = measurement_names.index(measurement)
|
|
data = fits.getdata(job)
|
|
assert data.shape[0] == expected_dims
|
|
assert data.shape[1] == expected_dims
|
|
filtered_data = np.zeros_like(data)
|
|
above_thresh_indexes = gaussian(data, 3) > thresholds[idx]
|
|
filtered_data[above_thresh_indexes] = 1
|
|
ratio_above_thresh = np.count_nonzero(above_thresh_indexes) / data.shape[0] / data.shape[1]
|
|
|
|
if ratio_above_thresh > ratio_above_thresh_max:
|
|
print(f"Exceeded ratio_above_thresh_max, possible data corruption in file: {job}")
|
|
new_name = job.split(".fits")[0] + "_e.fits"
|
|
os.rename(job, new_name)
|
|
|
|
# plt.figure(f"Data {measurement}")
|
|
# plt.imshow(data, cmap='jet')
|
|
# plt.figure(f"Circle")
|
|
# plt.imshow(filtered_data, cmap='jet')
|
|
# plt.show()
|
|
|
|
continue
|
|
|
|
hough_radii = np.arange(387, 394, 1)
|
|
hough_res = hough_circle(filtered_data, hough_radii)
|
|
|
|
vals, cxs, cys, rads = hough_circle_peaks(hough_res, hough_radii, threshold=circle_hough_thresh, total_num_peaks=10)
|
|
found_circle = None
|
|
for v, cx, cy, rad in zip(vals, cxs, cys, rads):
|
|
xskew = abs(639 - cx)
|
|
yskew = abs(639 - cy)
|
|
if xskew < max_center_skew and yskew < max_center_skew:
|
|
found_circle = (cx, cy, rad)
|
|
break
|
|
|
|
if not found_circle:
|
|
print(f"Could not find valid solar disc in file: {job}")
|
|
|
|
plt.figure(f"Data {measurement}")
|
|
if found_circle: cv.circle(data, (int(found_circle[0]),int(found_circle[1])), int(found_circle[2]), float(np.max(np.max(data))), 1)
|
|
plt.imshow(data, cmap='jet')
|
|
plt.figure(f"Circle")
|
|
if found_circle: cv.circle(filtered_data, (int(found_circle[0]),int(found_circle[1])), int(found_circle[2]), 0.5, 1)
|
|
plt.imshow(filtered_data, cmap='jet')
|
|
plt.show()
|
|
|
|
new_name = job.split(".fits")[0] + "_e.fits"
|
|
os.rename(job, new_name)
|
|
continue
|
|
|
|
# We have validated this file, rename it appropriately
|
|
new_name = job.split(".fits")[0] + "_f.fits"
|
|
os.rename(job, new_name)
|
|
|
|
|
|
except KeyboardInterrupt:
|
|
return
|
|
except Exception as e:
|
|
print(f"Error on file: {job} - {e}")
|
|
traceback.print_exception(e)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
stored_fits_dirs = [r"..\Data\goes18\l2\data\suvi-l2-ci094\2024",
|
|
r"..\Data\goes18\l2\data\suvi-l2-ci131\2024",
|
|
r"..\Data\goes18\l2\data\suvi-l2-ci171\2024",
|
|
r"..\Data\goes18\l2\data\suvi-l2-ci195\2024",
|
|
r"..\Data\goes18\l2\data\suvi-l2-ci284\2024",
|
|
r"..\Data\goes18\l2\data\suvi-l2-ci304\2024",]
|
|
# stored_fits_dirs = [r"..\Data\goes16\l2\data"]
|
|
# stored_fits_dirs = [r"Z:\NOAA GOES Data\fits_test_2024"]
|
|
|
|
reprocess_errors = False
|
|
nworkers = 20
|
|
|
|
|
|
regex_filename = r"dr_suvi-l2-ci\d{3}_g(16|18)_s\S*\.fits$"
|
|
|
|
|
|
work_queue = Queue(maxsize = nworkers)
|
|
workers = []
|
|
for i in range(nworkers):
|
|
p = Process(target = filter_fits, args = (work_queue,), daemon=True)
|
|
p.start()
|
|
workers.append(p)
|
|
|
|
lowpriority()
|
|
files_to_process = []
|
|
try:
|
|
for stored_fits_dir in stored_fits_dirs:
|
|
filename_tester = re.compile(regex_filename)
|
|
|
|
files_sorted_by_timestamp = defaultdict(list)
|
|
found_files = 0
|
|
print(f"Searching for FITS files in: {stored_fits_dir}")
|
|
for root, dirs, files in tqdm.tqdm(os.walk(stored_fits_dir), desc="Searching"):
|
|
for f in files:
|
|
if filename_tester.match(f):
|
|
file_parts = f.split("_")
|
|
file_name_end = file_parts[-1].split(".")[0]
|
|
if file_name_end == "f":
|
|
# if file_parts[-2] == "f":
|
|
# new_file_name = "_".join(file_parts[:-2]) + "_f.fits"
|
|
# print("bad file rename, fixing")
|
|
# os.rename(os.path.join(root,f), os.path.join(root,new_file_name))
|
|
|
|
# Already filtered and known good
|
|
# else:
|
|
# new_file_name = "_".join(file_parts[:-1]) + ".fits"
|
|
# print("Removing filter check")
|
|
# os.rename(os.path.join(root,f), os.path.join(root,new_file_name))
|
|
continue
|
|
elif file_name_end == "e":
|
|
if reprocess_errors:
|
|
new_file_name = "_".join(file_parts[:-1]) + ".fits"
|
|
os.rename(os.path.join(root,f), os.path.join(root,new_file_name))
|
|
files_to_process.append(os.path.abspath(os.path.join(root, new_file_name)))
|
|
else:
|
|
continue
|
|
elif file_name_end == "v1-0-2":
|
|
files_to_process.append(os.path.abspath(os.path.join(root,f)))
|
|
else:
|
|
print(f"Error - Unexpeted FITS file name: {f}")
|
|
for ftp in tqdm.tqdm(files_to_process, desc="Filtering files"):
|
|
work_queue.put(ftp)
|
|
except KeyboardInterrupt:
|
|
print("Finishing current jobs and exiting")
|
|
|
|
|
|
for _ in range(nworkers):
|
|
try:
|
|
work_queue.put(None, timeout=10.0)
|
|
except:
|
|
break
|
|
|
|
for w in workers:
|
|
w.join()
|