Merge pull request #11 from greglucas/site-rotation

Add site_rotation to add geographic ability to waveforms
greglucas · Dec 6, 2023 · d83d60a · d83d60a
2 parents 904fbd9 + c8054b5
commit d83d60a
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 0 deletions.
diff --git a/bezpy/mt/io.py b/bezpy/mt/io.py
@@ -109,6 +109,17 @@ def read_xml(fname):
     site.elevation = convert_float(get_text(loc, "Elevation"))
     site.declination = convert_float(get_text(loc, "Declination"))
 
+    # Orientation of the channels to geographic north (angle in degrees)
+    site.channel_orientation = {}
+    channel_mapping = {"Hx": "Bx", "Hy": "By", "Hz": "Bz", "Ex": "Ex", "Ey": "Ey"}
+    site_layout = root.find("SiteLayout")
+    for input_output in site_layout:
+        # InputChannels or OutputChannels subelements
+        for channel in input_output:
+            # Individual channels within the subelement
+            channel_name_mapped = channel_mapping[channel.attrib["name"]]
+            site.channel_orientation[channel_name_mapped] = convert_float(channel.attrib["orientation"])
+
     site.start_time = convert_datetime(get_text(xml_site, "Start"))
     site.end_time = convert_datetime(get_text(xml_site, "End"))
 

diff --git a/bezpy/mt/site.py b/bezpy/mt/site.py
@@ -230,10 +230,14 @@ def download_waveforms(self):
                                                                  self.start_time,
                                                                  self.end_time,
                                                                  network_code=self.network_code)
+        self._rotate_waveforms()
 
     def load_waveforms(self, directory="./"):
         """Load the waveform data that has already been downloaded."""
         self.waveforms = pd.read_hdf(directory + self.name + ".hdf", "waveforms")
+        if "Bx" not in self.waveforms.columns:
+            # Only rotate the waveforms if they haven't been rotated/don't exist yet
+            self._rotate_waveforms()
 
     def save_waveforms(self, directory="./"):
         """Save the waveform data to the specified file location."""
@@ -243,6 +247,24 @@ def save_waveforms(self, directory="./"):
                              ".download_waveforms()")
         self.waveforms.to_hdf(directory + self.name + ".hdf", "waveforms")
 
+    def _rotate_waveforms(self):
+        """Rotate the waveforms to the geographic coordinate system.
+
+        This is added as an extra method to add the components to the in-memory
+        objects, rather than saving the transformed values to disk.
+        """
+        # Magnetic field components rotated by the channel orientation
+        self.waveforms["Bx"] = (self.waveforms["BN"] * np.cos(np.deg2rad(self.channel_orientation["Bx"]))
+                                + self.waveforms["BE"] * np.cos(np.deg2rad(self.channel_orientation["By"])))
+        self.waveforms["By"] = (self.waveforms["BN"] * np.sin(np.deg2rad(self.channel_orientation["Bx"]))
+                                + self.waveforms["BE"] * np.sin(np.deg2rad(self.channel_orientation["By"])))
+
+        # Electric field components rotated by the channel orientation
+        self.waveforms["Ex"] = (self.waveforms["EN"] * np.cos(np.deg2rad(self.channel_orientation["Ex"]))
+                                + self.waveforms["EE"] * np.cos(np.deg2rad(self.channel_orientation["Ey"])))
+        self.waveforms["Ey"] = (self.waveforms["EN"] * np.sin(np.deg2rad(self.channel_orientation["Ex"]))
+                                + self.waveforms["EE"] * np.sin(np.deg2rad(self.channel_orientation["Ey"])))
+
 
 class Site1d(Site):
     """MT class for a 1D Fernberg profile"""