jeremy/noaa-goes-visualization
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

filter now processes images in time order and allows filtering by time

Browse source
This commit is contained in:
Jeremy Karst 2024-06-04 17:17:26 -04:00
parent bd225d1c83
commit 4471e1c503
1 changed files with 40 additions and 32 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

72
filter_FITS.py
View file

@ -3,6 +3,8 @@ from collections import defaultdict
from multiprocessing import Queue, Process from multiprocessing import Queue, Process
import re import re
import traceback import traceback
import datetime
import calendar
import tqdm import tqdm
import numpy as np import numpy as np
@ -42,7 +44,7 @@ def lowpriority():
measurement_names = ["094", "131", "171", "195", "284", "304"] measurement_names = ["094", "131", "171", "195", "284", "304"]
thresholds = [0.050, 0.10 , 1.00, 1.40, 0.90, 2.50] thresholds = [0.050, 0.10 , 1.00, 1.40, 0.90, 2.50]
max_center_skew = 7 max_center_skew = 7
circle_hough_thresh = 0.7 circle_hough_thresh = 0.75
expected_dims = 1280 expected_dims = 1280
ratio_above_thresh_max = 0.5 ratio_above_thresh_max = 0.5
@ -76,7 +78,7 @@ def filter_fits(work_queue):
continue continue
hough_radii = np.arange(387, 394, 1) hough_radii = np.arange(383, 393, 1)
hough_res = hough_circle(filtered_data, hough_radii) 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) 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) found_circle = (cx, cy, rad)
break 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: if not found_circle:
print(f"Could not find valid solar disc in file: {job}") 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" new_name = job.split(".fits")[0] + "_e.fits"
os.rename(job, new_name) os.rename(job, new_name)
continue continue
@ -116,23 +118,25 @@ def filter_fits(work_queue):
if __name__ == "__main__": if __name__ == "__main__":
stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci094", stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci094\2023",
r"..\Data\goes16\l2\data\suvi-l2-ci131", r"..\Data\goes16\l2\data\suvi-l2-ci131\2023",
r"..\Data\goes16\l2\data\suvi-l2-ci171", r"..\Data\goes16\l2\data\suvi-l2-ci171\2023",
r"..\Data\goes16\l2\data\suvi-l2-ci195", r"..\Data\goes16\l2\data\suvi-l2-ci195\2023",
r"..\Data\goes16\l2\data\suvi-l2-ci284", r"..\Data\goes16\l2\data\suvi-l2-ci284\2023",
r"..\Data\goes16\l2\data\suvi-l2-ci304", r"..\Data\goes16\l2\data\suvi-l2-ci304\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci094", r"..\Data\goes18\l2\data\suvi-l2-ci094\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci131", r"..\Data\goes18\l2\data\suvi-l2-ci131\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci171", r"..\Data\goes18\l2\data\suvi-l2-ci171\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci195", r"..\Data\goes18\l2\data\suvi-l2-ci195\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci284", r"..\Data\goes18\l2\data\suvi-l2-ci284\2023",
r"..\Data\goes18\l2\data\suvi-l2-ci304",] r"..\Data\goes18\l2\data\suvi-l2-ci304\2023",]
# stored_fits_dirs = [r"..\Data\goes16\l2\data"] # stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci131\2023\06\24"]
# stored_fits_dirs = [r"Z:\NOAA GOES Data\fits_test_2024"]
reprocess_errors = False starttime = calendar.timegm(datetime.datetime(2023, 1, 1, tzinfo=datetime.timezone.utc).timetuple())
nworkers = 16 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$" 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) workers.append(p)
lowpriority() lowpriority()
files_to_process = [] files_by_timestamp = defaultdict(list)
try: try:
for stored_fits_dir in stored_fits_dirs: for stored_fits_dir in stored_fits_dirs:
filename_tester = re.compile(regex_filename) filename_tester = re.compile(regex_filename)
files_sorted_by_timestamp = defaultdict(list)
found_files = 0 found_files = 0
print(f"Searching for FITS files in: {stored_fits_dir}") 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 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): if filename_tester.match(f):
file_parts = f.split("_") file_parts = f.split("_")
file_name_end = file_parts[-1].split(".")[0] 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": if file_name_end == "f":
continue # Already filtered from a previous run continue # Already filtered from a previous run
elif file_name_end == "e": elif file_name_end == "e":
if reprocess_errors: if reprocess_errors:
new_file_name = "_".join(file_parts[:-1]) + ".fits" new_file_name = "_".join(file_parts[:-1]) + ".fits"
os.rename(os.path.join(root,f), os.path.join(root,new_file_name)) 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: else:
continue continue
elif file_name_end == "v1-0-2": # This is the normal case for unprocessed data 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: else:
# print(f"Error - Unexpected FITS file name: {f}") # print(f"Error - Unexpected FITS file name: {f}")
pass pass
for ftp in tqdm.tqdm(files_to_process, desc="Filtering files"): sorted_times = sorted(list(files_by_timestamp.keys()))
work_queue.put(ftp) for st in tqdm.tqdm(sorted_times, desc="Filtering files"):
for f in files_by_timestamp[st]:
work_queue.put(f)
except KeyboardInterrupt: except KeyboardInterrupt:
print("Finishing current jobs and exiting") print("Finishing current jobs and exiting")