From 90d93b87a529747a7dfca6ec681cc5a34b955e07 Mon Sep 17 00:00:00 2001
From: Derisis13 <email.felhasznalasra@gmail.com>
Date: Sun, 6 Oct 2024 01:23:41 +0200
Subject: [PATCH] feat: detection on training sets, works with great accuracy

---
 create_matched_filter.py | 36 ++++++++++++++++++++++++++++++++++++
 detect_moonquakes.py     | 40 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 76 insertions(+)
 create mode 100644 create_matched_filter.py
 create mode 100644 detect_moonquakes.py

diff --git a/create_matched_filter.py b/create_matched_filter.py
new file mode 100644
index 0000000..a4bb0f6
--- /dev/null
+++ b/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()
diff --git a/detect_moonquakes.py b/detect_moonquakes.py
new file mode 100644
index 0000000..ebdb6ea
--- /dev/null
+++ b/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()