feat: detection on training sets, works with great accuracy

create_matched_filter.py

@@ -0,0 +1,36 @@
+"""Create the matched filter for correlational detection"""
+from datetime import datetime
+import numpy as np
+import matplotlib.pyplot as plt
+
+from bisect_moonquake import CAT_LUNAR, PREPROCESSED_LUNAR_DIR, from_mseed
+
+FILTER_SHAPE = [46376 // 3 + 1,] # I want it that way
+
+def create_matched_filter():
+    matched_filter = np.zeros(FILTER_SHAPE, 'O')
+    type_collection = CAT_LUNAR
+    for row in type_collection.iloc:
+        arrival_time = datetime.strptime(row['time_abs(%Y-%m-%dT%H:%M:%S.%f)'],'%Y-%m-%dT%H:%M:%S.%f')
+        sample_filename = row.filename + "_trimmed_7000_sec"
+        try:
+            st, _ = from_mseed(sample_filename, PREPROCESSED_LUNAR_DIR, arrival_time)
+        except FileNotFoundError:
+            # Because csv is faulty...
+            sample_filename = sample_filename.replace('HR00', 'HR02')
+            st, _ = from_mseed(sample_filename, PREPROCESSED_LUNAR_DIR, arrival_time)
+        st.traces[0].data *= 1 / st.traces[0].data.max()
+        st.traces[0].decimate(3)
+        try:
+            matched_filter += st.traces[0].data
+        except:
+            # print(st.traces[0].data.shape)
+            pass
+        # # Plot trace
+        # fig,ax = plt.subplots(1,1,figsize=(10,3))
+        # ax.plot(matched_filter)
+        # plt.show()
+    return matched_filter
+
+if __name__ == "__main__":
+    create_matched_filter()

detect_moonquakes.py

@@ -0,0 +1,40 @@
+import os
+from scipy.signal import correlate
+from bisect_moonquake import CAT_LUNAR, LUNAR_DATA_DIR, PREPROCESSED_LUNAR_DIR, from_mseed
+from datetime import datetime
+import create_matched_filter
+import matplotlib.pyplot as plt
+from pathlib import Path
+from tqdm import tqdm
+
+IMGDIR = "./images"
+def main():
+    matched_filter = create_matched_filter.create_matched_filter()
+    outfolder = Path(IMGDIR + "/lunar/training/")
+    outfolder.mkdir(parents=True, exist_ok=True)
+    for row in tqdm(CAT_LUNAR.iloc):
+        arrival_time = datetime.strptime(row['time_abs(%Y-%m-%dT%H:%M:%S.%f)'],'%Y-%m-%dT%H:%M:%S.%f')
+        test_filename = row.filename
+        try:
+            st, arrival = from_mseed(test_filename, LUNAR_DATA_DIR, arrival_time)
+        except FileNotFoundError:
+            # Because csv is faulty...
+            test_filename = test_filename.replace('HR00', 'HR02')
+            st, arrival = from_mseed(test_filename, LUNAR_DATA_DIR, arrival_time)
+        st.traces[0].data *= 1 / st.traces[0].data.max()
+        st.traces[0].decimate(3)
+        likelyhood = correlate(st.traces[0].data, matched_filter, mode='same')
+        estimated_arrival = (likelyhood.argmax() - matched_filter.shape[0] / 2) * st.traces[0].stats.delta
+        print(arrival - estimated_arrival)
+        # Plot trace
+        outfile = outfolder.joinpath(test_filename + "_correlation.svg")
+        fig,ax = plt.subplots(1,1,figsize=(10,3))
+        ax.plot(st.traces[0].times(), likelyhood)
+        ax.axvline(x = arrival, color='red',label='Rel. Arrival')
+        ax.axvline(x = estimated_arrival, color='green',label='Est. Arrival')
+        # plt.show()
+        plt.savefig(outfile)
+
+
+if __name__ == "__main__":
+    main()