From 4471e1c5037465ed2729dcb5780e76df28113293 Mon Sep 17 00:00:00 2001 From: Jeremy Karst Date: Tue, 4 Jun 2024 17:17:26 -0400 Subject: [PATCH] filter now processes images in time order and allows filtering by time --- filter_FITS.py | 72 ++++++++++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 32 deletions(-) diff --git a/filter_FITS.py b/filter_FITS.py index 6c0735d..7b7eec2 100644 --- a/filter_FITS.py +++ b/filter_FITS.py @@ -3,6 +3,8 @@ from collections import defaultdict from multiprocessing import Queue, Process import re import traceback +import datetime +import calendar import tqdm import numpy as np @@ -42,7 +44,7 @@ def lowpriority(): 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 +circle_hough_thresh = 0.75 expected_dims = 1280 ratio_above_thresh_max = 0.5 @@ -76,7 +78,7 @@ def filter_fits(work_queue): continue - hough_radii = np.arange(387, 394, 1) + hough_radii = np.arange(383, 393, 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) @@ -88,17 +90,17 @@ def filter_fits(work_queue): found_circle = (cx, cy, rad) break + # 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: {found_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() + 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 @@ -116,23 +118,25 @@ def filter_fits(work_queue): if __name__ == "__main__": - stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci094", - r"..\Data\goes16\l2\data\suvi-l2-ci131", - r"..\Data\goes16\l2\data\suvi-l2-ci171", - r"..\Data\goes16\l2\data\suvi-l2-ci195", - r"..\Data\goes16\l2\data\suvi-l2-ci284", - r"..\Data\goes16\l2\data\suvi-l2-ci304", - r"..\Data\goes18\l2\data\suvi-l2-ci094", - r"..\Data\goes18\l2\data\suvi-l2-ci131", - r"..\Data\goes18\l2\data\suvi-l2-ci171", - r"..\Data\goes18\l2\data\suvi-l2-ci195", - r"..\Data\goes18\l2\data\suvi-l2-ci284", - r"..\Data\goes18\l2\data\suvi-l2-ci304",] - # stored_fits_dirs = [r"..\Data\goes16\l2\data"] - # stored_fits_dirs = [r"Z:\NOAA GOES Data\fits_test_2024"] + stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci094\2023", + r"..\Data\goes16\l2\data\suvi-l2-ci131\2023", + r"..\Data\goes16\l2\data\suvi-l2-ci171\2023", + r"..\Data\goes16\l2\data\suvi-l2-ci195\2023", + r"..\Data\goes16\l2\data\suvi-l2-ci284\2023", + r"..\Data\goes16\l2\data\suvi-l2-ci304\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci094\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci131\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci171\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci195\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci284\2023", + r"..\Data\goes18\l2\data\suvi-l2-ci304\2023",] + # stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci131\2023\06\24"] - reprocess_errors = False - nworkers = 16 + starttime = calendar.timegm(datetime.datetime(2023, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) + stoptime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) + + reprocess_errors = True + nworkers = 20 regex_filename = r"dr_suvi-l2-ci\d{3}_g(16|18)_s\S*\.fits$" @@ -146,12 +150,11 @@ if __name__ == "__main__": workers.append(p) lowpriority() - files_to_process = [] + files_by_timestamp = defaultdict(list) 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"): @@ -159,22 +162,27 @@ if __name__ == "__main__": if filename_tester.match(f): file_parts = f.split("_") file_name_end = file_parts[-1].split(".")[0] + measure_end_time = int(datetime.datetime.strptime(file_parts[4][1:16] + " +0000", "%Y%m%dT%H%M%S %z").timestamp()) + if (measure_end_time >= starttime) and (measure_end_time < stoptime): + files_by_timestamp[measure_end_time].append(os.path.join(root,f)) if file_name_end == "f": continue # Already filtered from a previous run 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))) + files_by_timestamp[measure_end_time].append(os.path.abspath(os.path.join(root, new_file_name))) else: continue elif file_name_end == "v1-0-2": # This is the normal case for unprocessed data - files_to_process.append(os.path.abspath(os.path.join(root,f))) + files_by_timestamp[measure_end_time].append(os.path.abspath(os.path.join(root,f))) else: # print(f"Error - Unexpected FITS file name: {f}") pass - for ftp in tqdm.tqdm(files_to_process, desc="Filtering files"): - work_queue.put(ftp) + sorted_times = sorted(list(files_by_timestamp.keys())) + for st in tqdm.tqdm(sorted_times, desc="Filtering files"): + for f in files_by_timestamp[st]: + work_queue.put(f) except KeyboardInterrupt: print("Finishing current jobs and exiting")