output looks corrected but incorrectly rotated

generating_ASAGI_file.py

@@ -65,41 +65,42 @@ def writeNetcdf4SeisSol(sname, x, y, z, aName, aData):
     rootgrp.close()
 
 
-x = np.linspace(-50000, 50000, int((50000 + 50000) / 100) + 1)  # for testing if correct
-y = np.linspace(0, 10000, int(10000 / 100))
-z = np.linspace(-50000, 50000, int((50000 + 50000) / 100))
+if __name__ == '__main__':
+    x = np.linspace(-50000, 50000, int((50000 + 50000) / 100) + 1)  # for testing if correct
+    y = np.linspace(0, 10000, int(10000 / 100))
+    z = np.linspace(-50000, 50000, int((50000 + 50000) / 100))
 
-# Generate a Gaussian shaped slip distribution for illustration purposes
-xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')
-Rx = np.zeros(xg.shape)
-Rz = np.zeros(xg.shape)
+    # Generate a Gaussian shaped slip distribution for illustration purposes
+    xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')
+    Rx = np.zeros(xg.shape)
+    Rz = np.zeros(xg.shape)
 
-loc_temp = xg < -9800
-Rx[loc_temp] = (-xg[loc_temp] - 9800) / 10000
-loc_temp = xg > 1100
-Rx[loc_temp] = (xg[loc_temp] - 1100) / 10000
+    loc_temp = xg < -9800
+    Rx[loc_temp] = (-xg[loc_temp] - 9800) / 10000
+    loc_temp = xg > 1100
+    Rx[loc_temp] = (xg[loc_temp] - 1100) / 10000
 
-loc_temp = zg < -8000
-Rz[loc_temp] = (-zg[loc_temp] - 8000) / 1000
-loc_temp = zg > -2300
-Rz[loc_temp] = (zg[loc_temp] + 2300) / 10000
+    loc_temp = zg < -8000
+    Rz[loc_temp] = (-zg[loc_temp] - 8000) / 1000
+    loc_temp = zg > -2300
+    Rz[loc_temp] = (zg[loc_temp] + 2300) / 10000
 
-Rt = np.minimum(1, np.sqrt(np.power(Rx, 2) + np.power(Rz, 2)))
-s_xy0 = 30000000.0 * (1.0 - Rt)
+    Rt = np.minimum(1, np.sqrt(np.power(Rx, 2) + np.power(Rz, 2)))
+    s_xy0 = 30000000.0 * (1.0 - Rt)
 
-radius = np.sqrt(np.power(xg + 6000, 2) + np.power(zg + 6000, 2))
-pi = 4 * math.atan(1.0)
+    radius = np.sqrt(np.power(xg + 6000, 2) + np.power(zg + 6000, 2))
+    pi = 4 * math.atan(1.0)
 
-s_xy1 = np.zeros(xg.shape)
-loc_temp = radius <= 550
-s_xy1[loc_temp] = 3150000
-loc_temp = 1575000. * (1.0 + np.cos(pi * (radius - 550.0) / 250.0))
+    s_xy1 = np.zeros(xg.shape)
+    loc_temp = radius <= 550
+    s_xy1[loc_temp] = 3150000
+    loc_temp = 1575000. * (1.0 + np.cos(pi * (radius - 550.0) / 250.0))
 
-sheer_stress = s_xy0 + s_xy1
+    sheer_stress = s_xy0 + s_xy1
 
-prefix = "test"
-ldataName = ["sheer_stress"]
-lgridded_myData = [sheer_stress]
+    prefix = "test"
+    ldataName = ["sheer_stress"]
+    lgridded_myData = [sheer_stress]
 
-writeNetcdf4Paraview(prefix, x, y, z, ldataName, lgridded_myData)
-writeNetcdf4SeisSol(prefix, x, y, z, ldataName, lgridded_myData)
+    writeNetcdf4Paraview(prefix, x, y, z, ldataName, lgridded_myData)
+    writeNetcdf4SeisSol(prefix, x, y, z, ldataName, lgridded_myData)

get_point.py

@@ -3,19 +3,24 @@
 import typer
 from scipy.spatial import KDTree
 from typing_extensions import Annotated
+from generating_ASAGI_file import writeNetcdf4SeisSol, writeNetcdf4Paraview
 
 from meshgen import get_cartesian, apply_centering_points, apply_rotation_points
 
 
 def main(
         meta_file: Annotated[str, typer.Argument(help="Path for the meta file generated by meshgen script")],
-        lat: Annotated[str, typer.Argument(help="Latitude to convert to point on the mesh")],
-        long: Annotated[str, typer.Argument(help="Longitude to convert to point on the mesh")],
+        lat: Annotated[str, typer.Option(help="Latitude to convert to point on the mesh")] = 37.341206616740266,
+        long: Annotated[str, typer.Option(help="Longitude to convert to point on the mesh")] = -121.88075569799896,
         depth: Annotated[
             float, typer.Option(help="How deep should the point be from the sea level (in meters)")] = -1000.0,
         # from sea level in meters
         round_to_closest_point: Annotated[
-            bool, typer.Option(help="Should the point be rounded to the nearest point on the fault line")] = True
+            bool, typer.Option(help="Should the point be rounded to the nearest point on the fault line")] = True,
+        point_field_resolution: Annotated[
+            float, typer.Option(help="Resolution for netcdf files (in m)")] = 100,
+        output_distance_from_faults: Annotated[
+            str, typer.Option(help="Generate netcdf file for distance from fault points")] = None
 ):
     with h5py.File(meta_file, "r") as f:
         closest_point = np.array([long, lat, 0])
@@ -31,6 +36,28 @@ def main(
         cart = get_cartesian(lat_deg=closest_point[1], lon_deg=closest_point[0], alt=depth)
         cart = apply_rotation_points(cart, rotation_matrix)
         cart = apply_centering_points(cart, center)
+
+        if output_distance_from_faults is not None:
+            if f.get("bounding_box") is None:
+                print("Could not find bounding box in meta file, create bounding box with mesh")
+                exit()
+
+            print("generating KDTree")
+            fault_tree = KDTree(f["fault_points"][:])
+
+            bounding_box = f.get("bounding_box")[:]
+            min_coords = np.min(bounding_box, axis=0)
+            max_coords = np.max(bounding_box, axis=0)
+
+            x = np.linspace(min_coords[0], max_coords[0], int((max_coords[0] - min_coords[0]) / point_field_resolution) + 1)
+            y = np.linspace(min_coords[1], max_coords[1], int((max_coords[0] - min_coords[0]) / point_field_resolution)+ 2)
+            z = np.linspace(min_coords[2], max_coords[2], int((max_coords[0] - min_coords[0]) / point_field_resolution)+ 3)
+
+            xg, yg, zg = np.meshgrid(x, y, z, indexing='ij')
+            distances, index = fault_tree.query(np.stack((xg.flatten(), yg.flatten(), zg.flatten())).T, k=[1])
+            distances = distances.squeeze().reshape(xg.shape)
+
+            writeNetcdf4Paraview(output_distance_from_faults, x, y, z, ["fault_distance"], [distances])
         print(f"Location of the point is {cart}")