Skip to content

Commit

Permalink
removing trailing semicolons
Browse files Browse the repository at this point in the history
  • Loading branch information
kmaterna committed Dec 3, 2023
1 parent d760128 commit b2283be
Show file tree
Hide file tree
Showing 10 changed files with 537 additions and 526 deletions.
184 changes: 97 additions & 87 deletions src/InSAR_Profiles/grdtrack_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,130 +12,140 @@
class GMTProfile:
"""The internal object for a single GMT GRDTRACK profile."""
def __init__(self, center_lon, center_lat, azimuth, pt_lons, pt_lats, pt_dist, pt_azimuth, pt_value, offset=0):
self.center_lon = center_lon; self.center_lat = center_lat; self.azimuth = azimuth;
self.pt_lons = pt_lons; self.pt_lats = pt_lats; self.pt_dist = pt_dist;
self.pt_azimuth = pt_azimuth; self.pt_value = pt_value; self.offset = offset;
self.center_lon = center_lon # float
self.center_lat = center_lat # float
self.azimuth = azimuth # float
self.pt_lons = pt_lons # 1d array
self.pt_lats = pt_lats # 1d array
self.pt_dist = pt_dist # 1d array
self.pt_azimuth = pt_azimuth # 1d array
self.pt_value = pt_value # 1d array
self.offset = offset # float. essentially a type of metadata associated with each profile.

def set_offset(self, offset):
self.offset = offset;
self.offset = offset

def is_surface_breaking(self, critical_distance=0.25, critical_slope=0.30):
"""Is the slope at the center region of the profile sharp enough to be considered surface-breaking?"""
idx = np.where(np.abs(self.pt_dist) < critical_distance);
central_region = self.pt_value[idx];
derivative = np.diff(central_region);
idx = np.where(np.abs(self.pt_dist) < critical_distance)
central_region = self.pt_value[idx]
derivative = np.diff(central_region)
if sum(np.isnan(derivative)) == len(derivative):
return 0
if np.nanmax(np.abs(derivative)) > critical_slope:
return 1;
return 1
else:
return 0;
return 0

def get_offsets_surface_breaking_simple(self, medfilt_window=101, neg_sample_idx=-25, pos_sample_idx=25):
"""The methods section for determining slip amount when surface does break."""
negative_side = self.pt_value[np.where(self.pt_dist < 0)];
positive_side = self.pt_value[np.where(self.pt_dist > 0)];
negative_filt = scipy.signal.medfilt(negative_side, medfilt_window);
positive_filt = scipy.signal.medfilt(positive_side, medfilt_window);
offset = negative_filt[neg_sample_idx] - positive_filt[pos_sample_idx];
return offset;
"""The methods for determining slip amount when surface does break."""
negative_side = self.pt_value[np.where(self.pt_dist < 0)]
positive_side = self.pt_value[np.where(self.pt_dist > 0)]
negative_filt = scipy.signal.medfilt(negative_side, medfilt_window)
positive_filt = scipy.signal.medfilt(positive_side, medfilt_window)
offset = negative_filt[neg_sample_idx] - positive_filt[pos_sample_idx]
return offset

def get_offsets_subsurface_profile_simple(self, medfilt_window=101, critical_distance=0.22):
"""The methods section for determining slip amount when surface doesn't break."""
filtered = scipy.signal.medfilt(self.pt_value, medfilt_window); # FILTER THE PROFILE
filtered = general_utils.detrend_signal(self.pt_dist, filtered); # THEN WE DETREND THE SIGNAL
filtered_zoomed = filtered[np.where(np.abs(self.pt_dist) <= critical_distance)];
offset = np.nanmax(filtered_zoomed) - np.nanmin(filtered_zoomed); # gives a reasonable fit.
return offset;
"""The methods for determining slip amount when surface doesn't break."""
filtered = scipy.signal.medfilt(self.pt_value, medfilt_window) # FILTER THE PROFILE
filtered = general_utils.detrend_signal(self.pt_dist, filtered) # THEN WE DE-TREND THE SIGNAL
filtered_zoomed = filtered[np.where(np.abs(self.pt_dist) <= critical_distance)]
if sum(np.isnan(filtered_zoomed)) == len(filtered_zoomed):
return 0
offset = np.nanmax(filtered_zoomed) - np.nanmin(filtered_zoomed) # gives a reasonable fit.
return offset


def read_gmt_profiles(infile):
"""
Read a set of profiles automatically created from GMT GRDTRACK. File-IO function.
Returns a list of GMTProfile objects.
"""
profiles = [];
pt_lons, pt_lats, pt_dist, pt_azimuth, pt_value = [], [], [], [], [];
profiles = []
pt_lons, pt_lats, pt_dist, pt_azimuth, pt_value = [], [], [], [], []
with open(infile, 'r') as ifile:
print("Reading file %s " % infile);
print("Reading file %s " % infile)
for line in ifile:
if '>' in line:
temp = line.split();
center_lon = float(temp[6].split('/')[0]);
center_lat = float(temp[6].split('/')[1]);
azimuth = float(temp[7].split('=')[1]);
temp = line.split()
center_lon = float(temp[6].split('/')[0])
center_lat = float(temp[6].split('/')[1])
azimuth = float(temp[7].split('=')[1])
if len(pt_lons) > 0:
new_profile = GMTProfile(center_lon=center_lon, center_lat=center_lat, azimuth=azimuth,
pt_lons=np.array(pt_lons), pt_lats=np.array(pt_lats),
pt_dist=np.array(pt_dist), pt_azimuth=np.array(pt_azimuth),
pt_value=np.array(pt_value));
profiles.append(new_profile);
pt_lons, pt_lats, pt_dist, pt_azimuth, pt_value = [], [], [], [], [];
pt_value=np.array(pt_value))
profiles.append(new_profile)
pt_lons, pt_lats, pt_dist, pt_azimuth, pt_value = [], [], [], [], []
else:
temp = line.split();
pt_lons.append(float(temp[0]));
pt_lats.append(float(temp[1]));
pt_dist.append(float(temp[2]));
pt_azimuth.append(float(temp[3]));
pt_value.append(float(temp[4]));
print(len(profiles));
return profiles;
temp = line.split()
pt_lons.append(float(temp[0]))
pt_lats.append(float(temp[1]))
pt_dist.append(float(temp[2]))
pt_azimuth.append(float(temp[3]))
pt_value.append(float(temp[4]))
return profiles


def get_nearest_profile(profile_list, target_pt):
"""Return the profile whose center point is closest to a target point, such as a creepmeter location."""
distance_list = [];
"""Return the profile whose center point is closest to a target point, such as a creep-meter location."""
distance_list = []
for item in profile_list:
distance_list.append(haversine.distance((item.center_lat, item.center_lon), (target_pt[1], target_pt[0])));
return profile_list[np.argmin(distance_list)], np.argmin(distance_list);
distance_list.append(haversine.distance((item.center_lat, item.center_lon), (target_pt[1], target_pt[0])))
return profile_list[np.argmin(distance_list)], np.argmin(distance_list)


def visualize_offsets(profiles, outfile, vmin=0, vmax=30):
def visualize_offsets(profile_list, outfile, vmin=0, vmax=30):
"""Draw a plot with the fault trace color-coded by offset."""
plt.figure(dpi=300);
cm = plt.cm.get_cmap('hot');
colors = [x.offset for x in profiles]
sc = plt.scatter([x.center_lon for x in profiles], [x.center_lat for x in profiles], c=colors,
vmin=vmin, vmax=vmax, cmap=cm);
plt.colorbar(sc, label="LOS (mm)");
plt.savefig(outfile);
return;


def visualize_profile_displacement(profiles, outfile):
if len(profiles) > 0:
plt.figure(dpi=300);
for i, item in enumerate(profiles):
plt.plot(item.pt_dist, item.pt_value, '.');
plt.savefig(outfile);
return;


def draw_surface_trace(profiles, outfile, surface_breaking=()):
plt.figure();
plt.plot([x.center_lon for x in profiles], [x.center_lat for x in profiles], '.b', label='Profiles');
plt.figure(dpi=300)
cm = plt.cm.get_cmap('hot')
colors = [x.offset for x in profile_list]
sc = plt.scatter([x.center_lon for x in profile_list], [x.center_lat for x in profile_list], c=colors,
vmin=vmin, vmax=vmax, cmap=cm)
plt.colorbar(sc, label="LOS (mm)")
plt.savefig(outfile)
return


def visualize_profile_displacement(profile_list, outfile):
if len(profile_list) > 0:
plt.figure(dpi=300)
for i, item in enumerate(profile_list):
plt.plot(item.pt_dist, item.pt_value, '.')
plt.savefig(outfile)
return


def draw_surface_trace(profile_list, outfile, surface_breaking=()):
plt.figure()
plt.plot([x.center_lon for x in profile_list], [x.center_lat for x in profile_list], '.b', label='Profiles')
plt.plot([x.center_lon for x in surface_breaking], [x.center_lat for x in surface_breaking], '.r',
label='Surface breaking');
plt.legend();
plt.savefig(outfile);
return;
label='Surface breaking')
plt.legend()
plt.savefig(outfile)
return


def write_profile_offsets(profiles, outfile):
"""Write text files with calculated profile offsets."""
print("Writing file %s " % outfile);
ofile = open(outfile, 'w');
ofile.write("# center_lon center_lat azimuth offset\n");
for item in profiles:
ofile.write("%f %f %f %f\n" % (item.center_lon, item.center_lat, item.azimuth, item.offset) );
ofile.close();
return;
def write_profile_offsets(profile_list, outfile):
"""Write text files with user-calculated profile offsets."""
print("Writing file %s " % outfile)
ofile = open(outfile, 'w')
ofile.write("# center_lon center_lat azimuth offset\n")
for item in profile_list:
ofile.write("%f %f %f %f\n" % (item.center_lon, item.center_lat, item.azimuth, item.offset) )
ofile.close()
return


def read_profile_offsets(infile):
"""Read text files with calculated profile offsets."""
print("Reading file %s " % infile);
profiles = [];
[lons, lats, azimuths, offsets] = np.loadtxt(infile, skiprows=1, unpack=True);
"""Read text files with user-calculated profile offsets."""
print("Reading file %s " % infile)
profiles = []
[lons, lats, azimuths, offsets] = np.loadtxt(infile, skiprows=1, unpack=True)
for i in range(len(lons)):
new_profile = GMTProfile(center_lon=lons[i], center_lat=lats[i], azimuth=azimuths[i], offset=offsets[i],
pt_azimuth=(), pt_value=(), pt_lons=(), pt_lats=(), pt_dist=());
profiles.append(new_profile);
return profiles;
pt_azimuth=(), pt_value=(), pt_lons=(), pt_lats=(), pt_dist=())
profiles.append(new_profile)
return profiles
76 changes: 38 additions & 38 deletions src/Inversion/GF_element/GF_element.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,29 +29,29 @@ class GF_element:

def __init__(self, disp_points, param_name='', upper_bound=0, lower_bound=0, slip_penalty=0, units='',
fault_dict_list=(), points=()):
self.disp_points = disp_points;
self.param_name = param_name;
self.upper_bound = upper_bound;
self.lower_bound = lower_bound;
self.slip_penalty = slip_penalty;
self.units = units;
self.fault_dict_list = fault_dict_list;
self.points = points; # coordinates of surface trace of fault, if applicable
self.disp_points = disp_points
self.param_name = param_name
self.upper_bound = upper_bound
self.lower_bound = lower_bound
self.slip_penalty = slip_penalty
self.units = units
self.fault_dict_list = fault_dict_list
self.points = points # coordinates of surface trace of fault, if applicable

def set_param_name(self, param_name):
self.param_name = param_name;
self.param_name = param_name

def set_lower_bound(self, lower_bound):
self.lower_bound = lower_bound;
self.lower_bound = lower_bound

def set_upper_bound(self, upper_bound):
self.upper_bound = upper_bound;
self.upper_bound = upper_bound

def set_units(self, units_str):
self.units = units_str;
self.units = units_str

def set_slip_penalty(self, slip_penalty):
self.slip_penalty = slip_penalty;
self.slip_penalty = slip_penalty


def get_GF_rotation_element(obs_disp_points, ep, target_region=(-180, 180, -90, 90), rot_name=''):
Expand All @@ -64,22 +64,22 @@ def get_GF_rotation_element(obs_disp_points, ep, target_region=(-180, 180, -90,
:param rot_name: string, optional metadata for naming the rotation (ex: ocb_)
:returns: one GF_element with rotation displacements in x, y, and z directions
"""
rot_disp_p = [];
rot_disp_p = []
for obs_item in obs_disp_points:
coords = [obs_item.lon, obs_item.lat];
coords = [obs_item.lon, obs_item.lat]
if obs_item.is_within_bbox(target_region):
mult = 1;
mult = 1
else:
mult = 0;
response_to_rot = euler_pole.point_rotation_by_Euler_Pole(coords, ep);
mult = 0
response_to_rot = euler_pole.point_rotation_by_Euler_Pole(coords, ep)
response = Displacement_points(lon=obs_item.lon, lat=obs_item.lat, dE_obs=mult*response_to_rot[0],
dN_obs=mult*response_to_rot[1], dU_obs=mult*response_to_rot[2], Se_obs=0,
Sn_obs=0, Su_obs=0, meas_type=obs_item.meas_type, refframe=obs_item.refframe,
name=obs_item.name);
rot_disp_p.append(response);
name=obs_item.name)
rot_disp_p.append(response)
rotation_gf = GF_element(disp_points=rot_disp_p, param_name=rot_name, upper_bound=1, lower_bound=-1,
slip_penalty=0, units='deg/Ma');
return rotation_gf;
slip_penalty=0, units='deg/Ma')
return rotation_gf


def get_GF_rotation_elements(obs_disp_points, target_region=(-180, 180, -90, 90), rot_name=''):
Expand All @@ -95,12 +95,12 @@ def get_GF_rotation_elements(obs_disp_points, target_region=(-180, 180, -90, 90)
:returns: list of 3 GF_elements with rotation displacements in x, y, and z directions
"""
xresponse = get_GF_rotation_element(obs_disp_points, ep=[0, 0, 1], rot_name=rot_name + 'x_rot',
target_region=target_region); # X direction
target_region=target_region) # X direction
yresponse = get_GF_rotation_element(obs_disp_points, ep=[90, 0, 1], rot_name=rot_name + 'y_rot',
target_region=target_region); # Y direction
target_region=target_region) # Y direction
zresponse = get_GF_rotation_element(obs_disp_points, ep=[0, 89.99, 1], rot_name=rot_name + 'z_rot',
target_region=target_region); # Z direction
return [xresponse, yresponse, zresponse];
target_region=target_region) # Z direction
return [xresponse, yresponse, zresponse]


def get_GF_leveling_offset_element(obs_disp_points):
Expand All @@ -110,23 +110,23 @@ def get_GF_leveling_offset_element(obs_disp_points):
:param obs_disp_points: list of disp_point_objects
:returns: a list of 1 GF_element, or an empty list if there is no leveling in this dataset
"""
total_response_pts = [];
lev_count = 0;
total_response_pts = []
lev_count = 0
for item in obs_disp_points:
if item.meas_type == "leveling":
response = Displacement_points(lon=item.lon, lat=item.lat, dE_obs=0, dN_obs=0, dU_obs=1,
Se_obs=0, Sn_obs=0, Su_obs=0, meas_type=item.meas_type);
lev_count += 1;
Se_obs=0, Sn_obs=0, Su_obs=0, meas_type=item.meas_type)
lev_count += 1
else:
response = Displacement_points(lon=item.lon, lat=item.lat, dE_obs=0, dN_obs=0, dU_obs=0,
Se_obs=0, Sn_obs=0, Su_obs=0, meas_type=item.meas_type);
total_response_pts.append(response);
Se_obs=0, Sn_obs=0, Su_obs=0, meas_type=item.meas_type)
total_response_pts.append(response)
lev_offset_gf = GF_element(disp_points=total_response_pts, param_name='lev_offset',
upper_bound=1, lower_bound=-1, slip_penalty=0, units='m/yr');
upper_bound=1, lower_bound=-1, slip_penalty=0, units='m/yr')
if lev_count == 0:
return [];
return []
else:
return [lev_offset_gf];
return [lev_offset_gf]


def add_gfs(GFs_list):
Expand All @@ -136,7 +136,7 @@ def add_gfs(GFs_list):
:param GFs_list: list of GF_elements with the same modeled_disp_points lists
:returns: list of disp_points
"""
new_pts = GFs_list[0].disp_points;
new_pts = GFs_list[0].disp_points
for i in range(1, len(GFs_list)):
new_pts = dpo.utilities.add_disp_points(new_pts, GFs_list[i].disp_points);
return new_pts;
new_pts = dpo.utilities.add_disp_points(new_pts, GFs_list[i].disp_points)
return new_pts
14 changes: 7 additions & 7 deletions src/Inversion/GF_element/outputs.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,19 +11,19 @@ def visualize_GF_elements(GF_elements_list, outdir, exclude_list=()):
:param exclude_list: optional list of GF_element.param_name to exclude from visualizing
"""
if exclude_list == 'all':
return;
return
for GF_el in GF_elements_list:
if GF_el.param_name in exclude_list: # plot these elements separately, like individual CSZ patches
continue;
print(GF_el.param_name);
continue
print(GF_el.param_name)
if GF_el.param_name == "CSZ":
scale_arrow = (1.0, 0.010, "1 cm");
scale_arrow = (1.0, 0.010, "1 cm")
else:
scale_arrow = (1.0, 0.001, "1 mm");
scale_arrow = (1.0, 0.001, "1 mm")
library.plot_fault_slip.map_source_slip_distribution(GF_el.fault_dict_list, os.path.join(outdir, "gf_" +
GF_el.param_name + "_only.png"),
disp_points=GF_el.disp_points,
region=[-127, -119.7, 37.7, 43.3],
scale_arrow=scale_arrow,
v_labeling_interval=0.001);
return;
v_labeling_interval=0.001)
return
Loading

0 comments on commit b2283be

Please sign in to comment.