updated plots

Author: jhdorsett

.DS_Store

@@ -76,9 +76,9 @@ def load(filename):
     return(pickle.load(pkl_file))
 
 class Mesh:
-    def __init__(self):
+    def __init__(self,file):
 
-        data = np.loadtxt(params.gps_file,delimiter=',')
+        data = np.loadtxt(file,delimiter=',')
         data = data[::params.decimate,:]
         self.lon = data[:,0]
         self.lat = data[:,1]
@@ -224,7 +224,7 @@ def plot_vels(self, Ve=None, Vn=None, colormap='viridis',plot_centroid=False,res
 
         return(axes)
 
-    def plot(self,values=None, colormap='viridis', scale=None, cbar=True, lonlat=False,edges=False):
+    def plot(self,values=None, colormap='viridis', scale=None, cbar=True, lonlat=False,edges=False,borders=False):
 
         if values is None:
             values=[self.score]
@@ -248,6 +248,8 @@ def plot(self,values=None, colormap='viridis', scale=None, cbar=True, lonlat=Fal
 
         for i, (val, ax) in enumerate(zip(values, axes)):
             plt.subplot(1, n_vals, i + 1)
+
+
             if edges:
                 tripcolor = ax.tripcolor(triang, val, cmap=colormap, norm=scale,edgecolors='k')
             else:
@@ -265,7 +267,10 @@ def plot(self,values=None, colormap='viridis', scale=None, cbar=True, lonlat=Fal
                 else:
                     ax.plot([self.SegEnds[:, 0], self.SegEnds[:, 2]],
                             [self.SegEnds[:, 1], self.SegEnds[:, 3]], 'k-', linewidth=1.5)
-        
+            
+            if borders:
+                tools.plot_borders(ax,lonlat)
+
         if n_vals == 1:
                     return(axes[0])
         return axes
@@ -462,8 +467,7 @@ def Inversion(mesh):
         GExx = Gpoint.GExx
         GExy = Gpoint.GExy
         GEyy = Gpoint.GEyy
-
-
+    
     G = np.concatenate([GVe[mesh.ind,:],GVn[mesh.ind,:]])
 
     numobs = mesh.nodes.shape[0]*2

__pycache__/bForStrain.cpython-39.pyc

@@ -0,0 +1,16 @@
+Table S1: ITRF2014-PMM rotation pole cartesian components in DEG/My.
+----------------------------------------------------------------------------------------------------------
+   plate  omega_x   omega_y   omega_z
+          -----------DEG/My----------
+----------------------------------------------------------------------------------------------------------
+   ANTA   -0.0688   -0.0900    0.1874
+   ARAB    0.3205   -0.0378    0.4011
+   AUST    0.4194    0.3284    0.3375
+   EURA   -0.0235   -0.1476    0.2140
+   INDI    0.3205   -0.0014    0.4038
+   NAZC   -0.0925   -0.4290    0.4508
+   NOAM    0.0066   -0.1928   -0.0176
+   NUBI    0.0274   -0.1704    0.2037
+   PCFC   -0.1135    0.2907   -0.6025
+   SOAM   -0.0751   -0.0835   -0.0389
+   SOMA   -0.0336   -0.2206    0.2454

__pycache__/params.cpython-39.pyc

@@ -10,6 +10,7 @@ dependencies:
   - pyvista
   - scipy
   - pyproj
+  - geopandas
   - pip
   - pip:
     - triangle

__pycache__/tools.cpython-39.pyc

@@ -2,43 +2,40 @@
 
 gps_file = 'data/GNSS_velocities_alpides.txt'    
 creeping_file = []
-origin = [32.5, 35] 
+origin = [32.5, 35]
 lon_range = [-15,142]
-lat_range = [5, 60]
+lat_range = [5, 60] 
 nom_node_spacing = 400.
 
-## West US configuration
+### West US configuration
 
 gps_file = 'data/Zeng_vels_augmented_PA_JDF.txt'
 creeping_file = 'data/creeping_faults.txt'  
 origin = [-120,34] 
 lon_range = [-127, -96]
 lat_range = [26, 54]
-nom_node_spacing = 500.
+nom_node_spacing = 100.
 
-# mesh decimation factor (speeds up runtime, useful for debugging large scenarios)
-# set = 1 for no decimation
-decimate = 8
+### Meshing parameters
+# mesh decimation factor (speeds up runtime, useful for debugging large scenarios), set = 1 for no decimation
+decimate = 1
 
 # creeping patch length
 patchL = 15.
 
 #option to refine mesh in vicinity of GPS data
 # 0 = no refinement, 1-4 provides some
 refine_mesh = 0
-
-
-
 refine = 1
+
 nu = 0.25
 Gshear = 1
 
 betas = [40]
-
-
-uncertainty = True
+# Compute uncertainties? 
+uncertainty = False
 # number of realizations of strain rate for each beta value
-num = 10
+num = 25
 # relative weight on fitting creep rate data (creeping faults)
 Wcreep = 1
 # optional two-step minimization of strain rates below threshold value

bForStrain.py

@@ -4,6 +4,7 @@
 import params
 import pandas as pd
 import pyproj
+import geopandas as gpd
 
 ###############
 # Coordinate transforms defined using pyproj, called via wrapper functions
@@ -16,6 +17,23 @@ def llh2local(lons,lats):
 def local2llh(xs,ys): 
      return(np.array(local_to_llh.transform(xs,ys)).T)
 
+### Read in borders...
+# naturalearthdata.com data downloaded from https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_10m_admin_0_countries.geojson
+borders = gpd.read_file('data/ne_10m_admin_0_countries.geojson')
+borders_xy = borders.to_crs(local_crs)
+
+def plot_borders(ax,lonlat=False,width=0.5):
+     
+    xlim = ax.get_xlim()
+    ylim = ax.get_ylim()
+
+    if lonlat:
+        borders.boundary.plot(ax=ax, linewidth=width, color='black')
+    else:
+        borders_xy.boundary.plot(ax=ax, linewidth=width, color='black')
+        
+    ax.set_xlim(xlim)
+    ax.set_ylim(ylim)
 
 
 class build_G_pointforce: