Rewrote filter

ffmpeg_video.py

@@ -10,12 +10,12 @@ from sortedcontainers import SortedDict
 import tqdm
 from PIL import Image
 
-path_to_images = r"..\composite\goes16"
+path_to_images = r"..\composite\goes18"
-output_file = r"..\goes16_2024.mp4"
+output_file = r"..\goes18_2023.mp4"
-interp_file = r"..\goes16_2024_interp.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())
+starttime = calendar.timegm(datetime.datetime(2023, 1, 1, tzinfo=datetime.timezone.utc).timetuple())
-stoptime = calendar.timegm(datetime.datetime(2025, 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

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)
+            data = cv.resize(data, dsize=(expected_dims // 2, expected_dims // 2), interpolation=cv.INTER_LINEAR)
-            above_thresh_indexes = gaussian(data, 3) > thresholds[idx]
+            filtered_data = np.copy(data)
-            filtered_data[above_thresh_indexes] = 1
+            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)
+            if ratio_above_thresh < ratio_above_thresh_min:
-            hough_res = hough_circle(filtered_data, hough_radii)
+                print(f"Below ratio_above_thresh_min, possible bad data in file: {job}")
-
-            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}")
-
                 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):
+                plt.figure(figsize=[10.24, 7.68])
-                    cv.circle(filtered_data, (int(x),int(y)), int(r), 0.5, 1)
+                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
 
-            # We have validated this file, rename it appropriately
+            # Decide whether the data is valid based on solar disc placement and levels.
-            new_name = job.split(".fits")[0] + "_f.fits"
+            xavg = np.average(filtered_data, 0)
-            os.rename(job, new_name)
+            yavg = np.average(filtered_data, 1)
-            # print(f"Renamed: {job} -> {new_name}")
+            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')
+                
 
         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",
+    stored_fits_dirs = [r"..\Data\goes16\l2\data\suvi-l2-ci094\2024",
-                        r"..\Data\goes16\l2\data\suvi-l2-ci131\2023",
+                        r"..\Data\goes16\l2\data\suvi-l2-ci131\2024",
-                        r"..\Data\goes16\l2\data\suvi-l2-ci171\2023",
+                        r"..\Data\goes16\l2\data\suvi-l2-ci171\2024",
-                        r"..\Data\goes16\l2\data\suvi-l2-ci195\2023",
+                        r"..\Data\goes16\l2\data\suvi-l2-ci195\2024",
-                        r"..\Data\goes16\l2\data\suvi-l2-ci284\2023",
+                        r"..\Data\goes16\l2\data\suvi-l2-ci284\2024",
-                        r"..\Data\goes16\l2\data\suvi-l2-ci304\2023",
+                        r"..\Data\goes16\l2\data\suvi-l2-ci304\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci094\2023",
+                        r"..\Data\goes18\l2\data\suvi-l2-ci094\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci131\2023",
+                        r"..\Data\goes18\l2\data\suvi-l2-ci131\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci171\2023",
+                        r"..\Data\goes18\l2\data\suvi-l2-ci171\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci195\2023",
+                        r"..\Data\goes18\l2\data\suvi-l2-ci195\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci284\2023",
+                        r"..\Data\goes18\l2\data\suvi-l2-ci284\2024",
-                        r"..\Data\goes18\l2\data\suvi-l2-ci304\2023",]
+                        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())
+    starttime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple())
-    stoptime = 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

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())
+    starttime = calendar.timegm(datetime.datetime(2024, 1, 1, tzinfo=datetime.timezone.utc).timetuple())
-    stoptime = 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}")
+                                    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')}")
-                                    break
                         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

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)
@@ -247,8 +256,3 @@ if __name__ == "__main__":
 
     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}")