adding new reader functions for Noel's CSZ files

src/Inversion/GF_element/GF_element.py

@@ -27,17 +27,32 @@ class GF_element:
     :type points: np.array
     """
 
-    def __init__(self, disp_points, param_name, upper_bound, lower_bound, slip_penalty, units,
+    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.fault_dict_list = fault_dict_list;
         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;
+
+    def set_lower_bound(self, lower_bound):
+        self.lower_bound = lower_bound;
+
+    def set_upper_bound(self, upper_bound):
+        self.upper_bound = upper_bound;
+
+    def set_units(self, units_str):
+        self.units = units_str;
+
+    def set_slip_penalty(self, slip_penalty):
+        self.slip_penalty = slip_penalty;
+
 
 def get_GF_rotation_element(obs_disp_points, ep, target_region=(-180, 180, -90, 90), rot_name=''):
     """

src/Inversion/GF_element/readers_writers.py

@@ -1,8 +1,10 @@
-import numpy as np
-from Elastic_stresses_py.PyCoulomb import fault_slip_object as fso, disp_points_object as dpo
+from Elastic_stresses_py.PyCoulomb import fault_slip_object as fso
 from Elastic_stresses_py.PyCoulomb.disp_points_object.disp_points_object import Displacement_points
-from src.Inversion.GF_element import GF_element
-
+import Elastic_stresses_py.PyCoulomb.fault_slip_triangle as fst
+from Tectonic_Utils.geodesy import fault_vector_functions as fvf
+from .GF_element import GF_element
+import scipy.io
+import numpy as np
 
 def write_csz_dist_fault_patches(gf_elements, model_results_vector, outfile_gmt, outfile_txt):
     """Write out slip results for a distributed CSZ model into GMT format"""
@@ -22,11 +24,11 @@ def write_csz_dist_fault_patches(gf_elements, model_results_vector, outfile_gmt,
 
 def read_distributed_GF_static1d(gf_file, geom_file, latlonfile, latlonbox=(-127, -120, 38, 52), unit_slip=False):
     """
-    Read results of Fred's Static1D file (e.g., stat2C.outCascadia), getting partially into the GF_element object.
+    Read results of Fred's Static1D file (e.g., stat2C.outCascadia), getting into a minimal GF_element object.
     We also restrict the range of fault elements using a bounding box
     If unit_slip, we divide by the imposed slip rate to get a 1 cm/yr Green's Function.
     Returns a list of lists of disp_point objects, and a matching list of fault patches.
-    We get partially into the GF_element object; the rest of the way research-specific code in the Humboldt driver.
+    We get into a minimal GF_element object; the rest of the way research-specific code in the Humboldt driver.
     """
     fault_patches = fso.file_io.io_static1d.read_stat2C_geometry(geom_file);
     gps_disp_locs = fso.file_io.io_static1d.read_disp_points_from_static1d(latlonfile);
@@ -46,6 +48,7 @@ def read_distributed_GF_static1d(gf_file, geom_file, latlonfile, latlonbox=(-127
     counter = 0;
     norm_factor = 1;
     disp_points_all_patches, all_patches, given_slip = [], [], [];
+    GF_elements = [];
     for i in range(len(fault_patches)):
         if not fault_patches[i].is_within_bbox(latlonbox):  # can restrict to the southern patches if desired
             counter = counter + len(gps_disp_locs);
@@ -65,62 +68,77 @@ def read_distributed_GF_static1d(gf_file, geom_file, latlonfile, latlonbox=(-127
             counter = counter + 1;
             disp_points_one_patch.append(disp_point);
             if counter == len(gps_disp_locs):
-                break;
-        disp_points_all_patches.append(disp_points_one_patch);
+                break
         fault_slip_patch = fault_patches[i].change_fault_slip(fault_patches[i].slip * norm_factor);
-        all_patches.append(fault_slip_patch);
+
+        new_gf = GF_element(fault_dict_list=[fault_slip_patch], disp_points=disp_points_one_patch);
+        GF_elements.append(new_gf);
         given_slip.append(fault_patches[i].slip);  # in mm
 
-    return disp_points_all_patches, all_patches, given_slip;
+    return GF_elements, given_slip;
 
 
-def read_insar_greens_functions(gf_file, fault_patches, param_name='', lower_bound=0, upper_bound=0):
-    """
-    Read pre-computed green's functions in the format matching the write-function.
-    Currently, only works for triangles.
-
-    :param gf_file: string, filename
-    :param fault_patches: list of fault_slip_objects or fault_slip_triangles
-    :param param_name: string
-    :param lower_bound: float
-    :param upper_bound: float
+def extract_given_patch_helper(nodes, idx):
+    """ nodes = 1859 x 3.  idx = [a b c]."""
+    xs = [nodes[idx[0]-1][0], nodes[idx[1]-1][0], nodes[idx[2]-1][0], nodes[idx[0]-1][0]];
+    ys = [nodes[idx[0]-1][1], nodes[idx[1]-1][1], nodes[idx[2]-1][1], nodes[idx[0]-1][1]];
+    depths = [nodes[idx[0]-1][2], nodes[idx[1]-1][2], nodes[idx[2]-1][2], nodes[idx[0]-1][2]];
+    return xs, ys, depths;
+
+
+def read_GFs_matlab_CSZ(gf_file, patch_number):
     """
-    GF_elements = [];
-    gf_data_array = np.loadtxt(gf_file);
-    lons, lats = gf_data_array[:, 0], gf_data_array[:, 1];
-    model_disp_pts = [];
-    for tlon, tlat in zip(lons, lats):
-        mdp = Displacement_points(lon=tlon, lat=tlat, dE_obs=0, dN_obs=0, dU_obs=0, Se_obs=0, Sn_obs=0, Su_obs=0,
-                                  meas_type='insar');
-        model_disp_pts.append(mdp);
-
-    for i, tri in enumerate(fault_patches):  # right now, only the triangle version is written.
-        changed_slip = tri.change_fault_slip(rtlat=1, dipslip=0, tensile=0);  # triangle-specific
-        changed_slip = changed_slip.change_reference_loc();  # triangle-specific interface
-        index = i+2;  # moving to the correct column in the GF file, skipping lon and lat.
-        los_defo = gf_data_array[:, index];
-        model_disp_pts = dpo.utilities.set_east(model_disp_pts, los_defo);
-        GF_elements.append(GF_element.GF_element(disp_points=model_disp_pts, fault_dict_list=[changed_slip], units='m',
-                                                 param_name=param_name, lower_bound=lower_bound,
-                                                 upper_bound=upper_bound, slip_penalty=0));
-    return GF_elements;
-
-
-def write_insar_greens_functions(GF_elements, outfile):
+    Read the Green's functions for the CSZ calculated in Materna et al., 2019 by Noel Bartlow.
+    Returns a list of Green's Functions elements.
     """
-    Serialize a bunch of InSAR Green's Functions into written text file, in meters, with rows for each InSAR point:
-    For each observation point: lon, lat, dLOS1, dLOS2, dLOS3.... [n fault patches].
+    print("Reading file %s " % gf_file);
+    data_structure = scipy.io.loadmat(gf_file);
+    kern = data_structure['Kern'];  # 165 x 303 (E, N, U for each grid element);
+    fault_patches, nodes = read_matlab_CSZ_patches(gf_file);
+
+    # Discovering what's inside
+    for item in data_structure.keys():
+        print(item, np.shape(data_structure[item]));
+
+    num_gf_elements = len(fault_patches);
+    gf_elements = [];
+    for patch_num in range(num_gf_elements):
+        model_disp_points = [];
+        for i in range(len(data_structure['Lons'])):
+            new_item = Displacement_points(lon=data_structure['Lons'][i][0], lat=data_structure['Lats'][i][0],
+                                           dE_obs=kern[(3*i)][patch_number], dN_obs=kern[(3*i)+1][patch_number],
+                                           dU_obs=kern[(3*i)+2][patch_number]);
+            model_disp_points.append(new_item);  # Read model disp_points associated with one fault patch.
+        fault_patches[patch_number] = fault_patches[patch_number].change_fault_slip(1.0, 0, 0);
+
+        one_GF = GF_element(disp_points=model_disp_points, param_name=str(patch_number), units='meters',
+                            fault_dict_list=[fault_patches[patch_number]]);
+        gf_elements.append(one_GF);
+    return gf_elements;
+
 
-    :param GF_elements: list of GF_elements with InSAR displacements in disp_points.
-    :param outfile: string
+def read_matlab_CSZ_patches(gf_file):
     """
-    print("Writing file %s " % outfile);
-    ofile = open(outfile, 'w');
-    for i, pt in enumerate(GF_elements[0].disp_points):
-        ofile.write('%f %f ' % (pt.lon, pt.lat));
-        point_displacements = [GF_el.disp_points[i].dE_obs for GF_el in GF_elements];
-        for x in point_displacements:
-            ofile.write(str(x)+" ");
-        ofile.write("\n");
-    ofile.close();
-    return;
+    Read matlab file format in Materna et al., 2019.
+    Returns a list of triangular fault patches and a list of node points.
+    """
+    print("Reading file %s " % gf_file);
+    data_structure = scipy.io.loadmat(gf_file);
+    nodes = data_structure['nd_ll'];  # all the nodes for the entire CSZ, a big array from Canada to MTJ.
+    elements = data_structure['el'];  # elements
+
+    # Open all the fault patches
+    fault_patches = [];
+    for i, item in enumerate(elements):
+        xs, ys, depths = extract_given_patch_helper(nodes, item);
+        reflon, reflat = xs[0], ys[0];
+        v1x, v1y = fvf.latlon2xy_single(xs[0], ys[0], reflon, reflat);
+        v2x, v2y = fvf.latlon2xy_single(xs[1], ys[1], reflon, reflat);
+        v3x, v3y = fvf.latlon2xy_single(xs[2], ys[2], reflon, reflat);
+        new_ft = fst.fault_slip_triangle.TriangleFault(vertex1=[v1x*1000, v1y*1000, depths[0]*1000],
+                                                       vertex2=[v2x*1000, v2y*1000, depths[1]*1000],
+                                                       vertex3=[v3x*1000, v3y*1000, depths[2]*1000],
+                                                       lon=reflon, lat=reflat, depth=depths[0],
+                                                       rtlat_slip=0, dip_slip=0);
+        fault_patches.append(new_ft);
+    return fault_patches, nodes;

src/Inversion/GF_element/rw_insar_gfs.py

@@ -0,0 +1,56 @@
+import numpy as np
+from Elastic_stresses_py.PyCoulomb import disp_points_object as dpo
+from Elastic_stresses_py.PyCoulomb.disp_points_object.disp_points_object import Displacement_points
+from .GF_element import GF_element
+
+
+def read_insar_greens_functions(gf_file, fault_patches, param_name='', lower_bound=0, upper_bound=0):
+    """
+    Read pre-computed green's functions in the format matching the write-function.
+    Currently, only works for triangles.
+
+    :param gf_file: string, filename
+    :param fault_patches: list of fault_slip_objects or fault_slip_triangles
+    :param param_name: string
+    :param lower_bound: float
+    :param upper_bound: float
+    """
+    GF_elements = [];
+    gf_data_array = np.loadtxt(gf_file);
+    lons, lats = gf_data_array[:, 0], gf_data_array[:, 1];
+    model_disp_pts = [];
+    for tlon, tlat in zip(lons, lats):
+        mdp = Displacement_points(lon=tlon, lat=tlat, dE_obs=0, dN_obs=0, dU_obs=0, Se_obs=0, Sn_obs=0, Su_obs=0,
+                                  meas_type='insar');
+        model_disp_pts.append(mdp);
+
+    for i, tri in enumerate(fault_patches):  # right now, only the triangle version is written.
+        changed_slip = tri.change_fault_slip(rtlat=1, dipslip=0, tensile=0);  # triangle-specific
+        changed_slip = changed_slip.change_reference_loc();  # triangle-specific interface
+        index = i+2;  # moving to the correct column in the GF file, skipping lon and lat.
+        los_defo = gf_data_array[:, index];
+        model_disp_pts = dpo.utilities.set_east(model_disp_pts, los_defo);
+        GF_elements.append(GF_element(disp_points=model_disp_pts, fault_dict_list=[changed_slip], units='m',
+                                      param_name=param_name, lower_bound=lower_bound,
+                                      upper_bound=upper_bound, slip_penalty=0));
+    return GF_elements;
+
+
+def write_insar_greens_functions(GF_elements, outfile):
+    """
+    Serialize a bunch of InSAR Green's Functions into written text file, in meters, with rows for each InSAR point:
+    For each observation point: lon, lat, dLOS1, dLOS2, dLOS3.... [n fault patches].
+
+    :param GF_elements: list of GF_elements with InSAR displacements in disp_points.
+    :param outfile: string
+    """
+    print("Writing file %s " % outfile);
+    ofile = open(outfile, 'w');
+    for i, pt in enumerate(GF_elements[0].disp_points):
+        ofile.write('%f %f ' % (pt.lon, pt.lat));
+        point_displacements = [GF_el.disp_points[i].dE_obs for GF_el in GF_elements];
+        for x in point_displacements:
+            ofile.write(str(x)+" ");
+        ofile.write("\n");
+    ofile.close();
+    return;

src/bin/Projects/Calipatria/insar_gps_driver.py

@@ -11,6 +11,7 @@
 import Geodesy_Modeling.src.InSAR_1D_Object as InSAR_1D
 import Geodesy_Modeling.src.Inversion.GF_element.GF_element as GF_element
 import Geodesy_Modeling.src.Inversion.GF_element.readers_writers as GF_element_rw
+import src.Inversion.GF_element.rw_insar_gfs
 from Geodesy_Modeling.src.InSAR_1D_Object.class_model import InSAR_1D_Object
 import Tectonic_Utilities.Tectonic_Utils.seismo.moment_calculations as mo
 import numpy as np
@@ -63,8 +64,8 @@ def read_gf_elements_kalin(fault_file: str, gf_file: str):
     """Read a list of fault triangle elements and their associated GF's for use in inversion. """
     print("Reading pre-computed InSAR Green's functions.");
     fault_tris = fst.file_io.io_other.read_brawley_lohman_2005(fault_file);
-    GF_elements = GF_element_rw.read_insar_greens_functions(gf_file, fault_tris, param_name='kalin', lower_bound=-1,
-                                                            upper_bound=0);
+    GF_elements = src.Inversion.GF_element.rw_insar_gfs.read_insar_greens_functions(gf_file, fault_tris, param_name='kalin', lower_bound=-1,
+                                                                                    upper_bound=0);
     return GF_elements;
 
 
@@ -82,9 +83,9 @@ def write_misfit_report(exp_dict, obs_disp_pts, model_disp_pts):
 def compute_all_the_GFS(desc_pts: InSAR_1D_Object, asc_pts: InSAR_1D_Object):
     """ One-time function to compute and write ascending and descending green's functions (takes a few minutes)."""
     GF_elements_desc = compute_insar_gf_elements_kalin(exp_dict['fault_file'], desc_pts);
-    GF_element_rw.write_insar_greens_functions(GF_elements_desc, "desc_insar_gfs.txt");
+    src.Inversion.GF_element.rw_insar_gfs.write_insar_greens_functions(GF_elements_desc, "desc_insar_gfs.txt");
     GF_elements_asc = compute_insar_gf_elements_kalin(exp_dict['fault_file'], asc_pts);
-    GF_element_rw.write_insar_greens_functions(GF_elements_asc, "asc_insar_gfs.txt");
+    src.Inversion.GF_element.rw_insar_gfs.write_insar_greens_functions(GF_elements_asc, "asc_insar_gfs.txt");
     return;
 
 def combine_two_matching_lists_of_GF_elements(gf1, gf2):

src/bin/Projects/Humboldt/humboldt_inversion_driver.py

@@ -105,29 +105,26 @@ def read_hb_fault_gf_elements(exp_dict):
     for i in range(len(exp_dict["exp_faults"])):  # for each fault
         fault_name = exp_dict["exp_faults"][i];
         if fault_name == "CSZ_dist":  # Reading for distributed CSZ patches as unit slip.
-            one_patch_dps, csz_patches, maxslip = GF_element_rw.read_distributed_GF_static1d(exp_dict["inverse_dir"] +
-                                                                                             exp_dict["faults"]["CSZ"][
-                                                                                                 "GF"],
-                                                                                             exp_dict["inverse_dir"] +
-                                                                                             exp_dict["faults"]["CSZ"][
-                                                                                                 "geometry"],
-                                                                                             exp_dict["lonlatfile"],
-                                                                                             unit_slip=True,
-                                                                                             latlonbox=(
-                                                                                             -127, -120, 38, 44.5));
-            for gf_disp_points, patch, max0 in zip(one_patch_dps, csz_patches, maxslip):
-                lower_bound = exp_dict["faults"]["CSZ"]["slip_min"];  # default lower bound, probabaly zero
-                upper_bound = max0 * 130;  # upper bound about 40 mm/yr; from max_slip from geometry; units in cm
+            csz_GF_elements, maxslip = GF_element_rw.read_distributed_GF_static1d(exp_dict["inverse_dir"] +
+                                                                                  exp_dict["faults"]["CSZ"]["GF"],
+                                                                                  exp_dict["inverse_dir"] +
+                                                                                  exp_dict["faults"]["CSZ"]["geometry"],
+                                                                                  exp_dict["lonlatfile"],
+                                                                                  unit_slip=True,
+                                                                                  latlonbox=(-127, -120, 38, 44.5));
+            for gf_element, max0 in zip(csz_GF_elements, maxslip):  # experimental steps
+                gf_element.set_lower_bound(exp_dict["faults"]["CSZ"]["slip_min"]);  # default lower bound, probably zero
+                gf_element.set_upper_bound(max0 * 130);  # upper bound about 40 mm/yr; from geometry; units in cm
+                gf_element.set_param_name('CSZ_dist');
                 amount_of_slip_penalty = 1;
+                patch = gf_element.fault_dict_list[0];
                 if patch.depth > exp_dict["max_depth_csz_slip"]:
                     amount_of_slip_penalty = 100;  # optionally: force CSZ slip to be above a certain depth
                 if patch.depth < exp_dict["depth_of_forced_coupling"]:
-                    lower_bound = upper_bound * 0.90;  # optionally: force shallow CSZ to full coupling
-                one_gf_element = GF_element.GF_element(disp_points=gf_disp_points, param_name=fault_name,
-                                                       fault_dict_list=[patch], lower_bound=lower_bound,
-                                                       upper_bound=upper_bound, slip_penalty=amount_of_slip_penalty,
-                                                       units='cm/yr', points=());
-                gf_elements.append(one_gf_element);
+                    gf_element.set_lower_bound(gf_element.upper_bound * 0.90);  # optionally: force shallow CSZ coupled
+                gf_element.set_slip_penalty(amount_of_slip_penalty);
+                gf_element.set_units('cm/yr');
+                gf_elements.append(gf_element);
         else:  # Reading for LSF, MRF, other fault cases
             fault_gf = exp_dict["inverse_dir"] + exp_dict["faults"][fault_name]["GF"];
             fault_geom = exp_dict["inverse_dir"] + exp_dict["faults"][fault_name]["geometry"];