Gen seissol netcdf

generate_svg.py

@@ -30,8 +30,7 @@ def main(input_filename="outputs/BayModel1/merged_catmul_500m_a0.5.csv",
     proj = Proj(proj='utm', zone=get_utm_zone(bl_lon), ellps='WGS84', preserve_units=False)
     utm_bl_x, utm_bl_y = proj(bl_lon, bl_lat)
 
-
-    fig = plt.figure(frameon=False,figsize=(10,10*((utm_tr_x - utm_bl_x) / (utm_tr_y - utm_bl_y))))
+    fig = plt.figure(frameon=False, figsize=(5, 5 * ((utm_tr_x - utm_bl_x) / (utm_tr_y - utm_bl_y))))
     ax = fig.add_axes([0, 0, 1, 1])
 
     for line in data:
@@ -44,7 +43,7 @@ def main(input_filename="outputs/BayModel1/merged_catmul_500m_a0.5.csv",
             x, y = proj(lon_lat[0], lon_lat[1])
             utm_x.append(x)
             utm_y.append(y)
-            plt.plot(utm_x, utm_y, marker='o', linestyle='-')
+        plt.plot(utm_x, utm_y, marker='o', linestyle='-')
 
     plt.title('Projected UTM Coordinates')
     plt.xlabel('UTM Easting (m)')
@@ -56,7 +55,7 @@ def main(input_filename="outputs/BayModel1/merged_catmul_500m_a0.5.csv",
     # plt.gca().set_aspect((utm_tr_x - utm_bl_x) / (utm_tr_y - utm_bl_y))
 
     # Remove grid, legend, and borders
-    plt.axis('off')  # Turn off the axis
+    # plt.axis('off')  # Turn off the axis
 
     # Save the plot as an SVG file
     plt.savefig('utm_coordinates_plot.svg', format='svg')

get_point.py

@@ -32,6 +32,29 @@ def createNetcdf4ParaviewHandle(sname, x, y, z, aName):
     return rootgrp, vTd
 
 
+def createNetcdf4SeisSolHandle(sname, x, y, z, aName):
+    "create a netcdf file readable by paraview (but not by ASAGI)"
+    fname = sname + "_ASAGI.nc"
+    # print("writing " + fname)
+    # Creating the netcdf file
+    nx = x.shape[0]
+    ny = y.shape[0]
+    nz = z.shape[0]
+    rootgrp = Dataset(fname, "w", format="NETCDF4")
+    rootgrp.createDimension("u", nx)
+    rootgrp.createDimension("v", ny)
+    rootgrp.createDimension("w", nz)
+
+    vx = rootgrp.createVariable("u", "f4", ("u",))
+    vx[:] = x
+    vy = rootgrp.createVariable("v", "f4", ("v",))
+    vy[:] = y
+    vz = rootgrp.createVariable("w", "f4", ("w",))
+    vz[:] = z
+    vTd = rootgrp.createVariable("data", "f4", ("u", "v", "w"))
+    return rootgrp, vTd
+
+
 def main(
         meta_file: Annotated[str, typer.Argument(help="Path for the meta file generated by meshgen script")],
         lat: Annotated[str, typer.Option(help="Latitude to convert to point on the mesh")] = 37.341206616740266,
@@ -86,6 +109,7 @@ def main(
                 fault_tree = KDTree(f["fault_points"][:])
 
                 rootgrp, vTd = createNetcdf4ParaviewHandle(output_distance_from_faults, z, y, x, "fault_distance")
+                rootgrp_ss, vTd_ss = createNetcdf4SeisSolHandle(output_distance_from_faults, z, y, x, "fault_distance")
                 total_iterations = (len(x) // chunk_size + 1) * (len(y) // chunk_size + 1) * (len(z) // chunk_size + 1)
                 progress_bar = tqdm(total=total_iterations, desc="Generating fault_distance")
 
@@ -107,19 +131,22 @@ def main(
 
                             # Store the results in the appropriate slice of the distances array
                             vTd[k:k + chunk_size, j:j + chunk_size, i:i + chunk_size] = sub_distances
+                            vTd_ss[k:k + chunk_size, j:j + chunk_size, i:i + chunk_size] = sub_distances
 
                             # Update tqdm progress bar
                             progress_bar.update(1)
 
                 # Close tqdm progress bar
                 progress_bar.close()
                 rootgrp.close()
+                rootgrp_ss.close()
 
             if output_distance_from_topo is not None:
                 print("generating topography KDTree")
                 topo_tree = KDTree(f["topo_points"][:])
 
                 rootgrp, vTd = createNetcdf4ParaviewHandle(output_distance_from_topo, z, y, x, "topo_distance")
+                rootgrp_ss, vTd_ss = createNetcdf4SeisSolHandle(output_distance_from_topo, z, y, x, "topo_distance")
                 total_iterations = (len(x) // chunk_size + 1) * (len(y) // chunk_size + 1) * (len(z) // chunk_size + 1)
                 progress_bar = tqdm(total=total_iterations, desc="Generating topo_distance")
 
@@ -140,13 +167,15 @@ def main(
 
                             # Store the results in the appropriate slice of the distances array
                             vTd[k:k + chunk_size, j:j + chunk_size, i:i + chunk_size] = sub_distances
+                            vTd_ss[k:k + chunk_size, j:j + chunk_size, i:i + chunk_size] = sub_distances
 
                             # Update tqdm progress bar
                             progress_bar.update(1)
 
                 # Close tqdm progress bar
                 progress_bar.close()
                 rootgrp.close()
+                rootgrp_ss.close()
 
 
 if __name__ == "__main__":