From 4af3dc27ad44f17a9ab07a6f976cb157422b7a0d Mon Sep 17 00:00:00 2001 From: kmaterna Date: Mon, 17 Jun 2024 17:09:04 -0600 Subject: [PATCH] adding new experiment with synthetic creepmeter and LOS --- geodesy_modeling/Creepmeter/__init__.py | 0 .../Creepmeter/creepmeter_geometry.py | 30 +++++ .../Creepmeter/slip_1d_profiles.py | 41 ++++++ .../bin/Projects/Creep_Synthetics/driver.py | 118 ++++++++++++++++++ 4 files changed, 189 insertions(+) create mode 100644 geodesy_modeling/Creepmeter/__init__.py create mode 100644 geodesy_modeling/Creepmeter/creepmeter_geometry.py create mode 100644 geodesy_modeling/Creepmeter/slip_1d_profiles.py create mode 100755 geodesy_modeling/bin/Projects/Creep_Synthetics/driver.py diff --git a/geodesy_modeling/Creepmeter/__init__.py b/geodesy_modeling/Creepmeter/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/geodesy_modeling/Creepmeter/creepmeter_geometry.py b/geodesy_modeling/Creepmeter/creepmeter_geometry.py new file mode 100644 index 0000000..3ccf6df --- /dev/null +++ b/geodesy_modeling/Creepmeter/creepmeter_geometry.py @@ -0,0 +1,30 @@ +# Get the locations of points adjacent to the fault, i.e., a synthetic creep-meter + +import Tectonic_Utils.geodesy.haversine as haversine +import Tectonic_Utils.geodesy.fault_vector_functions as fault_vector_functions + + +def get_synthetic_creepmeter_points(lon0, lat0, strike, length_km, creepmeter_length): + """Calculate the two points that would make a synthetic creepmeter orthogonal to a particular fault. + The coordinates lon0 and lat0 are the back edge corner of the fault at the surface, not the center of the fault. + The total length is creepmeter_length (given in m)""" + + strike_vector = fault_vector_functions.get_strike_vector(strike) + west_vector = fault_vector_functions.get_strike_vector(strike - 90) + east_vector = fault_vector_functions.get_strike_vector(strike + 90) + + center_lon, center_lat = haversine.add_vector_to_coords(lon0, lat0, (length_km/2) * strike_vector[0], + (length_km/2) * strike_vector[1]) + west_coord = haversine.add_vector_to_coords(center_lon, center_lat, dx=(creepmeter_length/2000) * west_vector[0], + dy=(creepmeter_length/2000) * west_vector[1]) + east_coord = haversine.add_vector_to_coords(center_lon, center_lat, dx=(creepmeter_length/2000) * east_vector[0], + dy=(creepmeter_length/2000) * east_vector[1]) + return west_coord, east_coord + + +def write_cm_points(west_coord, east_coord, filename): + """Write the simple two points into a text file for espy calculations""" + with open(filename, 'w') as ofile: + ofile.write("%.9f %.9f \n" % (west_coord[0], west_coord[1])) + ofile.write("%.9f %.9f \n" % (east_coord[0], east_coord[1])) + return diff --git a/geodesy_modeling/Creepmeter/slip_1d_profiles.py b/geodesy_modeling/Creepmeter/slip_1d_profiles.py new file mode 100644 index 0000000..319449a --- /dev/null +++ b/geodesy_modeling/Creepmeter/slip_1d_profiles.py @@ -0,0 +1,41 @@ + +import numpy as np +import matplotlib.pyplot as plt + +""" +Assumes a 1d profile of slip in the following format: + +# Top_depth(km) bottom_depth(km) slip(mm). +2.50 3.50 0 +3.50 4.00 50 +4.00 5.00 0 +""" + + +def read_1d_profile(filename): + print("Reading filename %s" % filename) + top_depth, bottom_depth, slip_mm = np.loadtxt(filename, skiprows=1, unpack=True) + return top_depth, bottom_depth, slip_mm + + +def plot_1d_profile(txt_filename, plot_filename): + tops, bottoms, slips = read_1d_profile(txt_filename) + plt.figure(figsize=(4, 5), dpi=300) + x_vals, y_vals = [], [] + for x1, x2, slip in zip(tops, bottoms, slips): + x_vals.append(slip) + y_vals.append(x1) + x_vals.append(slip) + y_vals.append(x2) + plt.plot(x_vals, y_vals, linewidth=3, color='black') + plt.plot([0, np.max(x_vals)], [0, 0], linewidth=0.3, color='black') + plt.ylim([0, 4.3]) + plt.xlim([-1, 38]) + plt.xlabel('Slip (mm)', fontsize=14) + plt.ylabel('Depth (km)', fontsize=14) + plt.gca().invert_yaxis() + plt.gca().tick_params(axis='both', labelsize=14) + plt.tight_layout() + plt.savefig(plot_filename) + plt.close() + return diff --git a/geodesy_modeling/bin/Projects/Creep_Synthetics/driver.py b/geodesy_modeling/bin/Projects/Creep_Synthetics/driver.py new file mode 100755 index 0000000..83e0a74 --- /dev/null +++ b/geodesy_modeling/bin/Projects/Creep_Synthetics/driver.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python + +import os +import sys +from elastic_stresses_py.PyCoulomb import configure_calc, input_values, run_dc3d, output_manager +from geodesy_modeling.InSAR_2D_Object import model_enu_grids_into_los +from geodesy_modeling.Creepmeter import slip_1d_profiles, creepmeter_geometry + + +# 29 May 2024 +# Steps: +# Read 1D depth profile; plot it +# Build dislocation file, assuming a certain length and the creepmeter is in the middle of the dislocation +# Run Okada +# Predict displacements at creepmeter site +# Predict InSAR displacements, project into LOS. +# Plot InSAR displacements. + + +def configure_experiment(txt_filename): + _, filename = os.path.split(txt_filename) + exp_name = filename.split('.')[0] + os.makedirs("Outputs", exist_ok=True) + os.makedirs(os.path.join("Outputs", exp_name), exist_ok=True) + os.makedirs("espy_files", exist_ok=True) + return exp_name, os.path.join("Outputs", exp_name) + + +def build_okada_files(txt_filename, expname): + lon0, lat0, width_deg = -121.3, 36.7, 0.16 + strike, rake, dip, length_km = 321, 180, 89.9, 8 + tops, bottoms, slips = slip_1d_profiles.read_1d_profile(txt_filename) + minlon, maxlon, minlat, maxlat = lon0-width_deg, lon0+width_deg, lat0-width_deg, lat0+width_deg + espy_filename = os.path.join('espy_files', expname + '.intxt') + with open(espy_filename, 'w') as ofile: + print("Writing %s " % espy_filename) + ofile.write("# Top line\n") + ofile.write("General: 0.250 0.40 ") + ofile.write(str(minlon)+" "+str(maxlon)+" "+str(lon0)+" "+str(minlat)+" "+str(maxlat)+" "+str(lat0)+"\n") + for x1, x2, slip in zip(tops, bottoms, slips): + # Format: strike rake dip length_km width_km lon lat depth_km slip_m (opt: tensile_m) + width_km = x2 - x1 + ofile.write("Source_Patch: "+str(strike)+" "+str(rake)+" "+str(dip)+" "+str(length_km)+" ") + ofile.write(str(width_km)+' '+str(lon0)+' '+str(lat0)+' '+str(x1)+' '+str(slip/1000)+' 0\n') + + # Make the synthetic creepmeter locations + pt1, pt2 = creepmeter_geometry.get_synthetic_creepmeter_points(lon0, lat0, strike, length_km, creepmeter_length=10) + creepmeter_geometry.write_cm_points(pt1, pt2, filename=os.path.join('espy_files', 'target_pts.txt')) + return + + +def get_synthetic_creepmeter_magnitude(disp_points): + # Returns synthetic creepmeter displacements in mm + mag_v1 = disp_points[0].get_magnitude() + mag_v2 = disp_points[1].get_magnitude() + return 1000 * (mag_v1 + mag_v2) + + +def run_okada_calc(exp_name, outdir): + espy_filename = os.path.join('espy_files', exp_name + '.intxt') + params = configure_calc.Params(input_file=espy_filename, + disp_points_file=os.path.join("espy_files", "target_pts.txt"), + plot_stress=0, outdir=outdir) + [inputs, obs_disp_points, obs_strain_points] = input_values.read_inputs(params) + out_object = run_dc3d.do_stress_computation(params, inputs, obs_disp_points, obs_strain_points) + output_manager.produce_outputs(params, inputs, obs_disp_points, obs_strain_points, out_object) + synthetic_disp = get_synthetic_creepmeter_magnitude(out_object.model_disp_points) + print("SYNTHETIC DISP = ", synthetic_disp, " mm") + with open(os.path.join(outdir, 'synthetic_creepmeter.txt'), 'w') as ofile: + ofile.write("SYNTHETIC DISP = " + str(synthetic_disp) + " mm\n") + return out_object.model_disp_points + + +def project_to_los(outdir, disp_points=()): + los_params = { # REQUIRED PARAMETERS + 'east_grdfile': os.path.join(outdir, 'east.grd'), # with units of meters by default + 'north_grdfile': os.path.join(outdir, 'north.grd'), # with units of meters by default + 'up_grdfile': os.path.join(outdir, 'vert.grd'), # with units of meters by default + 'flight_direction': 190.0, + 'incidence_angle': 37.5, + 'wavelength_mm': 56, # with units of mm + 'look_direction': 'right', + 'plot_wrapped': True, # OPTIONAL PARAMETERS: + 'plot_unwrapped': True, + 'wrapped_plot_name': 'pred.jpg', + 'unwrapped_plot_name': 'unw_pred.jpg', + 'refidx': [0, 0], # either [float, float] for lon, lat, or [int, int] for row, col + 'plot_range': None, + 'plot_wrapped_annot': None, + 'plot_unwrapped_annot': None, + 'outdir': outdir, + 'label': 's1_descending_'} + model_enu_grids_into_los.do_synthetic_grid_LOS(los_params, disp_points=disp_points) + + los_params['flight_direction'] = 85 + los_params['incidence_angle'] = 52 + los_params['wavelength_mm'] = 240 + los_params['look_direction'] = 'left' + los_params['label'] = 'uavsar_eastward_' + model_enu_grids_into_los.do_synthetic_grid_LOS(los_params, disp_points=disp_points) + + los_params['flight_direction'] = 265 + los_params['incidence_angle'] = 52 + los_params['wavelength_mm'] = 240 + los_params['look_direction'] = 'left' + los_params['label'] = 'uavsar_westward_' + model_enu_grids_into_los.do_synthetic_grid_LOS(los_params, disp_points=disp_points) + return + + +if __name__ == "__main__": + given_filename = sys.argv[1] # the path to a 1d slip-distribution file + + exper_name, out_dir = configure_experiment(given_filename) + slip_1d_profiles.plot_1d_profile(given_filename, os.path.join(out_dir, 'slip_dist_' + exper_name + '.png')) + build_okada_files(given_filename, exper_name) + model_disp_points = run_okada_calc(exper_name, out_dir) + project_to_los(out_dir, model_disp_points)