Add Shapely Geometry Properties to Fault Objects

requirements.txt

@@ -5,9 +5,10 @@ pandas
 parse
 pygmt_helper @ git+https://github.com/ucgmsim/pygmt_helper
 qcore @ git+https://github.com/ucgmsim/qcore.git
+hypothesis[pandas, numpy]
 pytest
-hypothesis
 pytest-cov
+shapely
 scipy
 typer
-hypothesis[pandas, numpy]
+

source_modelling/sources.py

@@ -21,7 +21,7 @@
 
 import numpy as np
 import scipy as sp
-
+import shapely
 from qcore import coordinates, geo, grid
 
 _KM_TO_M = 1000
@@ -99,6 +99,11 @@ def centroid(self) -> np.ndarray:
         """np.ndarray: The centroid of the point source (which is just the point's coordinates)."""
         return self.coordinates
 
+    @property
+    def geometry(self) -> shapely.Point:
+        """shapely.Point: A shapely geometry for the point (projected onto the surface)."""
+        return shapely.Point(self.bounds)
+
     def fault_coordinates_to_wgs_depth_coordinates(
         self, fault_coordinates: np.ndarray
     ) -> np.ndarray:
@@ -272,7 +277,11 @@ def dip(self) -> float:
         """float: The dip angle of the fault."""
         return np.degrees(np.arcsin(np.abs(self.bottom_m - self.top_m) / self.width_m))
 
-    # TODO: change this to take width and dip rather than projected width.
+    @property
+    def geometry(self) -> shapely.Polygon:
+        """shapely.Polygon: A shapely geometry for the plane (projected onto the surface)."""
+        return shapely.Polygon(self.bounds)
+
     @staticmethod
     def from_centroid_strike_dip(
         centroid: np.ndarray,
@@ -501,6 +510,13 @@ def centroid(self) -> np.ndarray:
         """np.ndarray: The centre of the fault."""
         return self.fault_coordinates_to_wgs_depth_coordinates(np.array([1 / 2, 1 / 2]))
 
+    @property
+    def geometry(self) -> shapely.Polygon:
+        """shapely.Geometry: A shapely geometry for the fault (projected onto the surface)."""
+        return shapely.normalize(
+            shapely.union_all([plane.geometry for plane in self.planes])
+        )
+
     def wgs_depth_coordinates_to_fault_coordinates(
         self, global_coordinates: np.ndarray
     ) -> np.ndarray:

tests/test_sources.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import scipy as sp
+import shapely
 from hypothesis import assume, given, seed, settings
 from hypothesis import strategies as st
 from hypothesis.extra import numpy as nst
@@ -46,6 +47,10 @@ def test_point_construction(
     )
     assert np.allclose(point.coordinates, point_coordinates)
     assert point.length == point.width == length_m / 1000
+    assert np.allclose(
+        shapely.get_coordinates(point.geometry, include_z=True),
+        np.atleast_2d(point.bounds),
+    )
     assert np.isclose(point.width_m, point.width * 1000)
     assert np.allclose(point.centroid, point_coordinates)
 
@@ -148,6 +153,9 @@ def test_plane_construction(
     assert np.isclose(
         plane.width, plane.projected_width / np.cos(np.radians(plane.dip)), atol=1e-6
     )
+    assert np.allclose(
+        shapely.get_coordinates(plane.geometry, include_z=True)[:-1], plane.bounds
+    )
     assert np.isclose(plane.projected_width, projected_width, atol=1e-6)
     assert np.isclose(plane.strike, strike, atol=1e-6)
     assert np.isclose(plane.dip_dir, dip_dir, atol=1e-6)
@@ -248,6 +256,9 @@ def test_fault_construction(fault: Fault):
     assert np.isclose(fault.dip_dir, fault.planes[0].strike + 90)
     assert fault.corners.shape == (4 * len(fault.planes), 3)
     assert np.isclose(fault.area(), np.sum([plane.area for plane in fault.planes]))
+    assert fault.geometry.equals(
+        shapely.union_all([plane.geometry for plane in fault.planes])
+    )
     assert np.allclose(
         fault.wgs_depth_coordinates_to_fault_coordinates(fault.centroid),
         np.array([1 / 2, 1 / 2]),