From 756270cd1b7039cc2b863887e3213f775c7f465b Mon Sep 17 00:00:00 2001 From: Jeremy Karst Date: Fri, 5 Jul 2024 10:24:41 -0400 Subject: [PATCH] Rewrote filter --- ffmpeg_video.py | 10 +-- filter_FITS.py | 170 +++++++++++++++++++++++++++++++++--------------- merger_FITS.py | 13 ++-- puller_fits.py | 16 +++-- 4 files changed, 137 insertions(+), 72 deletions(-) diff --git a/ffmpeg_video.py b/ffmpeg_video.py index cff9720..1edf8bd 100644 --- a/ffmpeg_video.py +++ b/ffmpeg_video.py @@ -10,12 +10,12 @@ from sortedcontainers import SortedDict import tqdm from PIL import Image -path_to_images = r"..\composite\goes16" -output_file = r"..\goes16_2024.mp4" -interp_file = r"..\goes16_2024_interp.mp4" +path_to_images = r"..\composite\goes18" +output_file = r"..\goes18_2023.mp4" +interp_file = r"..\goes18_2023_interp.mp4" ffmpeg_path = r"..\ffmpeg.exe" -starttime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) -stoptime = calendar.timegm(datetime.datetime(2025, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) +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()) min_file_size = 350000 # Detect and remove corrupted files by filtering by file size max_file_size = 450000 encoding_crf = 16 diff --git a/filter_FITS.py b/filter_FITS.py index f9214a6..192078b 100644 --- a/filter_FITS.py +++ b/filter_FITS.py @@ -43,14 +43,22 @@ def lowpriority(): os.nice(19) measurement_names = ["094", "131", "171", "195", "284", "304"] -thresholds = [0.050, 0.10 , 1.00, 1.40, 0.95, 2.50] +thresholds = [0.050, 0.10 , 1.00, 1.40, 1.0, 2.50] max_center_skew = 7 -circle_hough_thresh = 0.75 expected_dims = 1280 +half_dims = expected_dims // 2 +valid_radii = [383, 394] +avg_radius = (valid_radii[0] + valid_radii[1]) // 2 +max_radius_error = 80 + ratio_above_thresh_max = 0.4 +ratio_above_thresh_min = 0.1 def filter_fits(work_queue): lowpriority() + idealcircle = np.zeros((expected_dims // 2, expected_dims // 2)) + cv.circle(idealcircle, (expected_dims // 4, expected_dims // 4), avg_radius // 2, 1, -1) + idealcircle_axis = np.average(idealcircle, 0) while True: try: job = work_queue.get() @@ -66,9 +74,11 @@ def filter_fits(work_queue): continue 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 + data = cv.resize(data, dsize=(expected_dims // 2, expected_dims // 2), interpolation=cv.INTER_LINEAR) + filtered_data = np.copy(data) + above_thresh_indexes = data > thresholds[idx] + filtered_data[above_thresh_indexes] = thresholds[idx] + filtered_data /= thresholds[idx] ratio_above_thresh = np.count_nonzero(above_thresh_indexes) / data.shape[0] / data.shape[1] if ratio_above_thresh > ratio_above_thresh_max: @@ -76,52 +86,104 @@ def filter_fits(work_queue): new_name = job.split(".fits")[0] + "_e.fits" os.rename(job, new_name) - plt.figure(f"Exceeded ratio_above_thresh_max", figsize=[10.24, 7.68]) + plt.figure(figsize=[10.24, 7.68]) + plt.title(f"Exceeded ratio_above_thresh_max [{ratio_above_thresh:0.2f}]") plt.imshow(filtered_data, cmap='jet') plt.savefig(job.split(".fits")[0] + "_e.jpg") plt.close('all') continue - hough_radii = np.arange(383, 394, 2) - 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=40) - 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 - - # print(f"Plotting file: {job}") - # plt.figure(f"Data {measurement}", figsize=[10.24, 7.68]) - # 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}", figsize=[10.24, 7.68]) - # 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}") - + if ratio_above_thresh < ratio_above_thresh_min: + print(f"Below ratio_above_thresh_min, possible bad data in file: {job}") new_name = job.split(".fits")[0] + "_e.fits" os.rename(job, new_name) - plt.figure(f"No Solar Disc Found", figsize=[10.24, 7.68]) - for x, y, r in zip(cxs, cys, rads): - cv.circle(filtered_data, (int(x),int(y)), int(r), 0.5, 1) + + plt.figure(figsize=[10.24, 7.68]) + plt.title(f"Below ratio_above_thresh_min [{ratio_above_thresh:0.2f}]") plt.imshow(filtered_data, cmap='jet') plt.savefig(job.split(".fits")[0] + "_e.jpg") plt.close('all') + + continue + + # Decide whether the data is valid based on solar disc placement and levels. + xavg = np.average(filtered_data, 0) + yavg = np.average(filtered_data, 1) + axis_indices = list(range(half_dims)) + + # Center Estimation + weighted_avg_center_x = np.average(axis_indices, 0, xavg) + weighted_avg_center_y = np.average(axis_indices, 0, yavg) + + good_center = True + if (half_dims // 2 - weighted_avg_center_x > max_center_skew) or (half_dims // 2 - weighted_avg_center_y > max_center_skew): + print(f"Could not find centered solar disc in file: {job}") + good_center = False + + # Radius Estimation + edge_thresh = 0.98 + high_edge_x = np.argmax(np.cumsum(xavg) > (np.sum(xavg) * edge_thresh)) + low_edge_x = len(xavg) - np.argmax(np.cumsum(np.flip(xavg)) > (np.sum(xavg) * edge_thresh)) + high_edge_y = np.argmax(np.cumsum(yavg) > (np.sum(yavg) * edge_thresh)) + low_edge_y = len(yavg) - np.argmax(np.cumsum(np.flip(yavg)) > (np.sum(yavg) * edge_thresh)) + + good_radius = True + calc_radius_x = (high_edge_x - low_edge_x) # We do not divide by two because our image is downsized by half in each dimension + calc_radius_y = (high_edge_y - low_edge_y) + calc_radius = (calc_radius_x + calc_radius_y) / 2.0 + if abs(calc_radius - avg_radius) > max_radius_error: + print(f"Could not find correct solar radius in file: {job}") + good_radius = False + + # Goodness of fit estimation vs perfect disc + gof_x = np.sum(np.abs(xavg - idealcircle_axis)) / half_dims + gof_y = np.sum(np.abs(yavg - idealcircle_axis)) / half_dims + + good_fit = True + if gof_x > 0.2 or gof_y > 0.2: + print(f"Could not find valid solar disc in file: {job}") + good_fit = False + + + if good_center and good_radius and good_fit: + # We have validated this file, rename it appropriately + new_name = job.split(".fits")[0] + "_f.fits" + os.rename(job, new_name) + else: + new_name = job.split(".fits")[0] + "_e.fits" + os.rename(job, new_name) + + # Plot Results + fig, axes = plt.subplots(2, 2, width_ratios=(0.2, 1), height_ratios=(0.2, 1), gridspec_kw={"hspace":0.0, "wspace": 0.0}, figsize = [8, 8]) + + axes[0][0].axis('off') + + axes[0][1].plot(axis_indices, xavg) + axes[0][1].plot(axis_indices, idealcircle_axis) + # axes[0][1].fill_between(axis_indices, xavg, idealcircle_axis, hatch="//", edgecolor="red", facecolor="none") + + axes[1][0].plot(yavg, axis_indices) + axes[1][0].plot(idealcircle_axis, axis_indices) + axes[1][0].invert_xaxis() + axes[1][0].invert_yaxis() + + axes[1][1].imshow(filtered_data, aspect='auto', vmin=0, vmax=1) + axes[1][1].vlines([low_edge_x, high_edge_x], 0, half_dims) + axes[1][1].hlines([low_edge_y, high_edge_y], 0, half_dims) + axes[1][1].scatter(half_dims // 2, half_dims // 2, marker="o", linewidths = 1, alpha = 0.75, edgecolors="black", facecolors = "none") + axes[1][1].scatter(weighted_avg_center_x, weighted_avg_center_y, marker = "x", alpha = 0.75, facecolors='red') + + # Set shared x axis between imshow and xdata plots + axes[0][1].xaxis.set_ticks_position("top") + axes[0][1].yaxis.set_ticks_position("right") + axes[0][1].sharex(axes[1][1]) + axes[1][0].sharey(axes[1][1]) + axes[1][1].yaxis.set_ticks_position("right") + + plt.savefig(job.split(".fits")[0] + "_e.jpg") + plt.close('all') - continue - - # We have validated this file, rename it appropriately - new_name = job.split(".fits")[0] + "_f.fits" - os.rename(job, new_name) - # print(f"Renamed: {job} -> {new_name}") except KeyboardInterrupt: return @@ -131,23 +193,23 @@ def filter_fits(work_queue): if __name__ == "__main__": - 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-ci094\2024", + r"..\Data\goes16\l2\data\suvi-l2-ci131\2024", + r"..\Data\goes16\l2\data\suvi-l2-ci171\2024", + r"..\Data\goes16\l2\data\suvi-l2-ci195\2024", + r"..\Data\goes16\l2\data\suvi-l2-ci284\2024", + r"..\Data\goes16\l2\data\suvi-l2-ci304\2024", + 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\goes18\l2\data\suvi-l2-ci284\2023\01\04"] + # stored_fits_dirs = [r"..\Data\goes18\l2\data\suvi-l2-ci284\2023\01\10"] - 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()) + starttime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) + stoptime = calendar.timegm(datetime.datetime(2025, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) reprocess_errors = False nworkers = 16 diff --git a/merger_FITS.py b/merger_FITS.py index 0c66553..833578d 100644 --- a/merger_FITS.py +++ b/merger_FITS.py @@ -435,13 +435,13 @@ image_names = ["094A", "131A", "171A", "195A", "284A", "304A"] if __name__ == "__main__": stored_fits_dirs = [r"..\Data\goes16\l2\data", r"..\Data\goes18\l2\data"] processed_images_dirs = [r"..\composite\goes16", r"..\composite\goes18"] - 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()) + starttime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) + stoptime = calendar.timegm(datetime.datetime(2025, 1, 1, tzinfo=datetime.timezone.utc).timetuple()) regex_filename = r"dr_suvi-l2-ci\d{3}_g(16|18)_s\S*\_f.fits$" - nworkers = 20 + nworkers = 16 max_time_gap = 3 - fill_missing_data = False + fill_missing_data = True file_prefixes = ["dr_suvi-l2-ci" + n[:-1] for n in image_names] @@ -533,7 +533,7 @@ if __name__ == "__main__": files_this_timestamp = sorted(files_this_timestamp) if (not len(files_this_timestamp) == 6): if fill_missing_data: - print(f"Invalid or incomplete sensor records for {timestamp} - {len(files_this_timestamp)}/6 filling from last good data.") + print(f"Invalid or incomplete sensor records for {timestamp} | {datetime.datetime.fromtimestamp(timestamp, datetime.UTC).strftime('%Y-%m-%d %H:%M:%S %Z')} - {len(files_this_timestamp)}/6 filling from last good data.") files_for_job = [] for i, prefix in enumerate(file_prefixes): found = False @@ -550,8 +550,7 @@ if __name__ == "__main__": if time_gap <= max_time_gap: files_for_job.append(last_good_files[i]) else: - print(f"Detected a gap of {time_gap} frames at {timestamp}") - break + print(f"Detected a gap of {time_gap} frames at {timestamp} | {datetime.datetime.fromtimestamp(timestamp, datetime.UTC).strftime('%Y-%m-%d %H:%M:%S %Z')}") if len(files_for_job) == 6: work_queue.put((files_for_job, timestamp, processed_images_dir, True)) else: # We have a complete file set, update the last_good_files diff --git a/puller_fits.py b/puller_fits.py index 78c99e4..6d3ba4c 100644 --- a/puller_fits.py +++ b/puller_fits.py @@ -10,6 +10,7 @@ from threading import Thread import queue import math from functools import partial +from queue import Empty import tqdm @@ -229,16 +230,24 @@ if __name__ == "__main__": break except KeyboardInterrupt: print("Saving file database and shutting down.") + except Empty: + print("Work Complete, Shutting down...") except Exception as e: print(f"Unhandled Exception during run: {e}") print("Shutting down") + with open(file_database_path, 'w') as f: + f.write(json.dumps(file_info_cache)) + for _ in range(nfetchworkers): try: query_work_queue.put(None, timeout=5.0) except: break + time.sleep(1) + print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") + for _ in range(ndownloadworkers): try: download_work_queue.put(None, timeout=5.0) @@ -246,9 +255,4 @@ if __name__ == "__main__": break for w in workers: - w.join(5.0) - - with open(file_database_path, 'w') as f: - f.write(json.dumps(file_info_cache)) - - print(f"Downloaded {fetched_image_count} | Already had {already_had_image_count} | Failed {failed_image_count}") \ No newline at end of file + w.join(5.0) \ No newline at end of file