diff --git a/README.md b/README.md index e495240..814a724 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,3 @@ -[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/Thomas-Ulrich/StressInversionNotebook/master) +[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/Thomas-Ulrich/StressInversionNotebook/HEAD?labpath=StressInversionNotebook.ipynb) # StressInversionNotebook This notebook was compiled by Thomas Ulrich and provides exercises on stress inversion. It is part of the Master's level "P7.2 Seismology" lecture at LMU Munich taught by Alice-Agnes Gabriel. diff --git a/StressInversionNotebook.ipynb b/StressInversionNotebook.ipynb index 78b0ca5..a23b551 100644 --- a/StressInversionNotebook.ipynb +++ b/StressInversionNotebook.ipynb @@ -30,9 +30,9 @@ "\n", "\n", "\n", - "The Wallace–Bott hypothesis, which relates **N** with **s**, allow relating $\\tau$ with **s**, and therefore with observed focal mechanisms. \n", + "The Wallace–Bott hypothesis, which relates **N** with **s**, allows relating $\\tau$ with **s**, and therefore with observed focal mechanisms. \n", " \n", - "In following derivation, we use the Einstein summation convention (an index variable appearing twice implies summation of that term over all the values of the index) and use the kronecker delta, defined by:\n", + "In the following derivation, we use the Einstein summation convention (an index variable appearing twice implies summation of that term over all the values of the index) and use the Kronecker delta, defined by:\n", "\n", "$\\delta_{mp}=1$ if m=p, 0 if m!=p\n", "\n", @@ -84,6 +84,8 @@ "import matplotlib.pyplot as plt\n", "from scipy.linalg import lstsq\n", "import seaborn as sns\n", + "import pandas as pd\n", + "\n", "sns.set()\n", "\n", "\n", @@ -92,12 +94,35 @@ ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "In the following cell, various focal mechanism datasets can be loaded depending on the `mydata` variable. Possible choices are `NZ` for the GeoNet database (see https://github.com/GeoNet/data/tree/master/moment-tensor), `Japan2000-2010` for a dataset extracted from the ISC catalog over a region englobing the location of the Tohoku-oki earthquake (see http://www.isc.ac.uk/iscbulletin/search/fmechanisms/#reviewed) in the period 2000-2010, and `Japan2012-2020` over the period 2012-2020. Other values for `mydata` loads the dataset from the stress inversion code `Stressinverse_1.1.2` over West Bohemia (Czech Republic).\n", + "def computeA(n1, n2, n3):\n", + " \"\"\"equation 8 in Vavrycuk et al. , GJI. 2014\"\"\"\n", + " A = np.zeros((3, 5))\n", + " n12 = n1**2\n", + " n22 = n2**2\n", + " n32 = n3**2\n", + " A[0, 0] = n1 * (n22 + 2 * n32)\n", + " A[0, 1] = n2 * (1 - 2 * n12)\n", + " A[0, 2] = n3 * (1 - 2 * n12)\n", + " A[0, 3] = n1 * (-n22 + n32)\n", + " A[0, 4] = -2 * n1 * n2 * n3\n", "\n", - "The `use_auxiliary_plane` variable can be set to `True` to check the consequence of focal mechanism wrongly set to the auxiliary plane rather than the fault plane." + " A[1, 0] = n2 * (-n12 + n32)\n", + " A[1, 1] = n1 * (1 - 2 * n22)\n", + " A[1, 2] = -2 * n1 * n2 * n3\n", + " A[1, 3] = n2 * (n12 + 2 * n32)\n", + " A[1, 4] = n3 * (1 - 2 * n22)\n", + "\n", + " A[2, 0] = n3 * (-2 * n12 - n22)\n", + " A[2, 1] = -2 * n1 * n2 * n3\n", + " A[2, 2] = n1 * (1 - 2 * n32)\n", + " A[2, 3] = n3 * (-n12 - 2 * n22)\n", + " A[2, 4] = n2 * (1 - 2 * n32)\n", + " return A" ] }, { @@ -106,109 +131,77 @@ "metadata": {}, "outputs": [], "source": [ - "# mydata = ''\n", - "# mydata = 'NZ'\n", - "mydata = \"Japan2000-2010\"\n", - "# mydata = 'Japan2012-2020'\n", + "def ComputeSlipVector(strike, dip, rake):\n", + " \"\"\"Compute slip vector from focal mechanisms\"\"\"\n", + " nobservations = strike.shape[0]\n", + " SlipVector = np.zeros((nobservations, 3))\n", "\n", - "# Just for testing\n", - "use_auxiliary_plane = False\n", + " # e.g. Stein and Wyssession (p218)\n", + " SlipVector[:, 0] = np.cos(rake[:]) * np.cos(strike[:]) + np.sin(rake[:]) * np.cos(\n", + " dip[:]\n", + " ) * np.sin(strike[:])\n", + " SlipVector[:, 1] = -np.cos(rake[:]) * np.sin(strike[:]) + np.sin(rake[:]) * np.cos(\n", + " dip[:]\n", + " ) * np.cos(strike[:])\n", + " SlipVector[:, 2] = np.sin(rake[:]) * np.sin(dip[:])\n", "\n", - "if mydata == \"NZ\":\n", - " # moment tensor solutions from Geonet are available for download here\n", - " # https://github.com/GeoNet/data/blob/master/moment-tensor/GeoNet_CMT_solutions.csv\n", - " # we download the raw data using:\n", - " # wget https://raw.githubusercontent.com/GeoNet/data/master/moment-tensor/GeoNet_CMT_solutions.csv\n", - " # HEADER\n", - " # 0 1 2 3 4 5 6 7 8 9 10 11 12 13\n", - " # PublicID,Date,Latitude,Longitude,strike1,dip1,rake1,strike2,dip2,rake2,ML,Mw,Mo,CD,NS,DC,Mxx,Mxy,Mxz,Myy,Myz,Mzz,VR,Tva,Tpl,Taz,Nva,Npl,Naz,Pva,Ppl,Paz\n", - "\n", - " data = np.genfromtxt(\"GeoNet_CMT_solutions.csv\", delimiter=\",\", skip_header=1)\n", - " lat = data[:, 2]\n", - " lon = data[:, 3]\n", - " strike = data[:, 4]\n", - " dip = data[:, 5]\n", - " rake = data[:, 6]\n", - " depth = data[:, 13]\n", - " # Auxiliary plane\n", - " strikeAP = data[:, 7]\n", - " dipAP = data[:, 8]\n", - " rakeAP = data[:, 9]\n", - "\n", - " # vizualize e.g. using https://www.openstreetmap.org. export->Manually select a different area\n", - "\n", - " # region around the Kaikoura earthquake\n", - " ids = np.where((lon > 172.7) & (lon < 174.9) & (lat > -42.7) & (lat < -40.5) & (depth < 20))[0]\n", + " return SlipVector.flatten()\n", "\n", - " # region around the Alpine Fault\n", - " # ids = np.where((lon > 166.6) & (lon < 170.6) & (lat > -44.6) & (lat < -43.3) & (depth < 20))[0]\n", "\n", - " # North of Northern island\n", - " # ids = np.where((lon > 176.1) & (lon < 178.9) & (lat > -40.3) & (lat < -36.8) & (depth > 15))[0]\n", + "def ComputeNormalVector(strike, dip, rake):\n", + " \"\"\"Compute fault normal from focal mechanisms\"\"\"\n", + " nobservations = strike.shape[0]\n", + " FaultNormal = np.zeros((nobservations, 3))\n", + " # e.g. Stein and Wyssession (p218)\n", + " FaultNormal[:, 0] = -np.sin(dip[:]) * np.sin(strike[:])\n", + " FaultNormal[:, 1] = -np.sin(dip[:]) * np.cos(strike[:])\n", + " FaultNormal[:, 2] = np.cos(dip[:])\n", "\n", - " # South of Southern island\n", - " # ids = np.where((lon > 163.1) & (lon < 169.5) & (lat > -49.1) & (lat < -44.9) & (depth > 15))[0]\n", + " return FaultNormal\n", "\n", - " strike = data[ids, 4]\n", - " dip = data[ids, 5]\n", - " rake = data[ids, 6]\n", "\n", - "elif mydata == \"Japan2000-2010\":\n", - " # Downloaded from http://www.isc.ac.uk/iscbulletin/search/fmechanisms/#reviewed\n", - " data = np.genfromtxt(\"Japan2000_2010.csv\", delimiter=\",\", skip_header=24)\n", - " lat = data[:, 4]\n", - " lon = data[:, 5]\n", - " depth = data[:, 6]\n", - " ids = np.where((lat < 38.2) & (lon > 141) & (depth > 20))\n", - " # ids = np.where(lon>140.5)\n", - " data = data[ids]\n", - " strike = data[:, 19]\n", - " dip = data[:, 20]\n", - " rake = data[:, 21]\n", - " # Auxiliary plane\n", - " strikeAP = data[:, 22]\n", - " dipAP = data[:, 23]\n", - " rakeAP = data[:, 24]\n", - "\n", - "elif mydata == \"Japan2012-2020\":\n", - " # Downloaded from http://www.isc.ac.uk/iscbulletin/search/fmechanisms/#reviewed\n", - " data = np.genfromtxt(\"Japan2012_2020.csv\", delimiter=\",\", skip_header=26)\n", - " lat = data[:, 4]\n", - " lon = data[:, 5]\n", - " depth = data[:, 6]\n", - " ids = np.where((lat < 38.2) & (lon > 141) & (depth > 20))\n", - " data = data[ids]\n", - " strike = data[:, 19]\n", - " dip = data[:, 20]\n", - " rake = data[:, 21]\n", - " # Auxiliary plane\n", - " strikeAP = data[:, 22]\n", - " dipAP = data[:, 23]\n", - " rakeAP = data[:, 24]\n", - "else:\n", - " # default file from Stressinverse_1.1.2\n", - " print(\"using default dataset\")\n", - " strike, dip, rake = np.loadtxt(\"West_Bohemia_mechanisms.dat\", skiprows=2, unpack=True)\n", + "def AssembleAmatrix(FaultNormal):\n", + " nobservations = FaultNormal.shape[0]\n", + " # Assemble A matrix\n", + " A = np.zeros((3 * nobservations, 5))\n", + " for i in range(nobservations):\n", + " A[3 * i : 3 * i + 3, :] = computeA(*FaultNormal[i, :])\n", + " return A\n", "\n", - "if use_auxiliary_plane:\n", - " # set for 30% of focal mechanism the wrong plane\n", - " nfocal = len(strike)\n", - " listn = np.random.choice(range(nfocal), size=int(0.3 * nfocal), replace=False)\n", - " strike[listn] = strikeAP[listn]\n", - " dip[listn] = dipAP[listn]\n", - " rake[listn] = rakeAP[listn]\n", "\n", - "# Remove nan if any\n", - "strike = strike[~np.isnan(strike)]\n", - "dip = dip[~np.isnan(dip)]\n", - "rake = rake[~np.isnan(rake)]\n", + "def computeStressTensor(A, S):\n", + " \"\"\"Solve the Problem At=s using the generalized inverse of A\n", + " A^-g = (A^t A)^-1 A^t (see e.g. Stein and Wyssession (p418, eq 17)\n", + " This the best solution in the least square sense\n", + " The textbook solution is:\n", + " Ag = A.T @ A\n", + " Ag_inv = np.linalg.pinv(Ag)\n", + " t = Ag_inv @ A.T @ S\n", + " but this squares the condition number of the matrix\n", + " Therefore we use a more robust function from scipy\n", + " \"\"\"\n", + " t, residual, rank, singular_values = lstsq(A, S)\n", + "\n", + " # t indices are 11, 12, 13, 22, 23\n", + " # stress tensor in ENU\n", + " # Tij = np.array([[t[0], t[1], t[2]], [t[1], t[3], t[4]], [t[2], t[4], -t[0]-t[3]]])\n", + " # stress tensor in ESD\n", + " Tij = np.array([[t[0], -t[1], -t[2]], [-t[1], t[3], t[4]], [-t[2], t[4], -t[0] - t[3]]])\n", + " return Tij\n", "\n", - "print(\"number of focal mechanism selected\", len(strike))\n", - "print(\"strike dip rake arrays:\\n\", strike, dip, rake)\n", "\n", - "strike0 = np.radians(strike)\n", - "dip0 = np.radians(dip)\n", - "rake0 = np.radians(rake)" + "def ComputeSortedEigenValuesAndVectors(Tij):\n", + " \"\"\"Compute eigen Values and eigen Vectors, to compute R ratio and SH_max\"\"\"\n", + "\n", + " # each column of eigenVectors is an eigenvector\n", + " eigenValues, eigenVectors = np.linalg.eig(Tij)\n", + " # important notes:\n", + " # the largest principal stress corresponds with the smallest seigeValue!\n", + " # the eigenValues should not be sorted in absolute value!\n", + " idx = eigenValues.argsort()\n", + " eigenValues = eigenValues[idx]\n", + " eigenVectors = eigenVectors[:, idx]\n", + " return eigenValues, eigenVectors" ] }, { @@ -217,30 +210,34 @@ "metadata": {}, "outputs": [], "source": [ - "def computeA(n1, n2, n3):\n", - " \"\"\"equation 8 in Vavrycuk et al. , GJI. 2014\"\"\"\n", - " A = np.zeros((3, 5))\n", - " n12 = n1 ** 2\n", - " n22 = n2 ** 2\n", - " n32 = n3 ** 2\n", - " A[0, 0] = n1 * (n22 + 2 * n32)\n", - " A[0, 1] = n2 * (1 - 2 * n12)\n", - " A[0, 2] = n3 * (1 - 2 * n12)\n", - " A[0, 3] = n1 * (-n22 + n32)\n", - " A[0, 4] = -2 * n1 * n2 * n3\n", + "def ComputeSHmax(eigenVectors):\n", + " # largest principal stress\n", + " uz = np.array([0.0, 0.0, 1.0])\n", + " idv = int(np.argmax(np.absolute(np.dot(eigenVectors.T, uz))))\n", + " # print(\"most vertical vector\", idv + 1)\n", "\n", - " A[1, 0] = n2 * (-n12 + n32)\n", - " A[1, 1] = n1 * (1 - 2 * n22)\n", - " A[1, 2] = -2 * n1 * n2 * n3\n", - " A[1, 3] = n2 * (n12 + 2 * n32)\n", - " A[1, 4] = n3 * (1 - 2 * n22)\n", + " # compute largest eigen value not the most vertical eigen vector\n", + " ih = 1 if idv == 0 else 0\n", "\n", - " A[2, 0] = n3 * (-2 * n12 - n22)\n", - " A[2, 1] = -2 * n1 * n2 * n3\n", - " A[2, 2] = n1 * (1 - 2 * n32)\n", - " A[2, 3] = n3 * (-n12 - 2 * n22)\n", - " A[2, 4] = n2 * (1 - 2 * n32)\n", - " return A" + " vh = eigenVectors[:, ih]\n", + " # SHmax: azimuth of the maximum horizontal compressive stress\n", + " SHmax = np.degrees(np.arctan2(vh[1], vh[0]))\n", + " return SHmax\n", + "\n", + "\n", + "def plotSHmaxR(eigenVectors, eigenValues):\n", + " # Compute SH_max and R (s2_ratio)\n", + " R = (eigenValues[:, 0] - eigenValues[:, 1]) / (eigenValues[:, 0] - eigenValues[:, 2])\n", + " plt.title(\"R distribution\", fontsize=16)\n", + " plt.hist(R)\n", + " plt.show()\n", + " nobservations = eigenValues.shape[0]\n", + " SHmax = np.zeros(nobservations)\n", + " for i in range(nobservations):\n", + " SHmax[i] = ComputeSHmax((eigenVectors[i, :, :]).T)\n", + " plt.title(r\"$SH_{max}$ distribution\", fontsize=16)\n", + " plt.hist(SHmax)\n", + " plt.show()" ] }, { @@ -377,7 +374,13 @@ "\n", " # --------------------------------------------------------------------------\n", " # legend\n", - " axSA.legend((sig1, sig2, sig3), (\"sigma 1\", \"sigma 2\", \"sigma 3\"), loc=\"lower right\", fontsize=14, numpoints=1)\n", + " axSA.legend(\n", + " (sig1, sig2, sig3),\n", + " (\"sigma 1\", \"sigma 2\", \"sigma 3\"),\n", + " loc=\"lower right\",\n", + " fontsize=14,\n", + " numpoints=1,\n", + " )\n", " if plot_file == \"\":\n", " plt.show()\n", " else:\n", @@ -385,79 +388,83 @@ " plt.close()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following cell, various focal mechanism datasets can be loaded depending on the `mydata` variable. Possible choices are `NZ` for the GeoNet database (see https://github.com/GeoNet/data/tree/master/moment-tensor), `Japan2000_2010` for a dataset extracted from the ISC catalog over a region englobing the location of the Tohoku-oki earthquake (see http://www.isc.ac.uk/iscbulletin/search/fmechanisms/#reviewed) in the period 2000-2010, and `Japan2012_2020` over the period 2012-2020. Other values for `mydata` loads the dataset from the stress inversion code `Stressinverse_1.1.2` over West Bohemia (Czech Republic).\n", + "\n", + "The `use_auxiliary_plane` variable can be set to `True` to check the consequence of focal mechanism wrongly set to the auxiliary plane rather than the fault plane." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "def ComputeSlipVector(strike, dip, rake):\n", - " \"\"\"Compute slip vector from focal mechanisms\"\"\"\n", - " nobservations = strike.shape[0]\n", - " SlipVector = np.zeros((nobservations, 3))\n", - "\n", - " # e.g. Stein and Wyssession (p218)\n", - " SlipVector[:, 0] = np.cos(rake[:]) * np.cos(strike[:]) + np.sin(rake[:]) * np.cos(dip[:]) * np.sin(strike[:])\n", - " SlipVector[:, 1] = -np.cos(rake[:]) * np.sin(strike[:]) + np.sin(rake[:]) * np.cos(dip[:]) * np.cos(strike[:])\n", - " SlipVector[:, 2] = np.sin(rake[:]) * np.sin(dip[:])\n", - "\n", - " return SlipVector.flatten()\n", - "\n", - "\n", - "def ComputeNormalVector(strike, dip, rake):\n", - " \"\"\"Compute fault normal from focal mechanisms\"\"\"\n", - " nobservations = strike.shape[0]\n", - " FaultNormal = np.zeros((nobservations, 3))\n", - " # e.g. Stein and Wyssession (p218)\n", - " FaultNormal[:, 0] = -np.sin(dip[:]) * np.sin(strike[:])\n", - " FaultNormal[:, 1] = -np.sin(dip[:]) * np.cos(strike[:])\n", - " FaultNormal[:, 2] = np.cos(dip[:])\n", - "\n", - " return FaultNormal\n", - "\n", - "\n", - "def AssembleAmatrix(FaultNormal):\n", - " nobservations = FaultNormal.shape[0]\n", - " # Assemble A matrix\n", - " A = np.zeros((3 * nobservations, 5))\n", - " for i in range(nobservations):\n", - " A[3 * i : 3 * i + 3, :] = computeA(*FaultNormal[i, :])\n", - " return A\n", - "\n", - "\n", - "def computeStressTensor(A, S):\n", - " \"\"\"Solve the Problem At=s using the generalized inverse of A\n", - " A^-g = (A^t A)^-1 A^t (see e.g. Stein and Wyssession (p418, eq 17)\n", - " This the best solution in the least square sense\n", - " The textbook solution is:\n", - " Ag = A.T @ A\n", - " Ag_inv = np.linalg.pinv(Ag)\n", - " t = Ag_inv @ A.T @ S\n", - " but this squares the condition number of the matrix\n", - " Therefore we use a more robust function from scipy\n", - " \"\"\"\n", - " t, residual, rank, singular_values = lstsq(A, S)\n", + "# mydata = ''\n", + "mydata = \"NZ\"\n", + "# mydata = \"Japan2000_2010\"\n", + "# mydata = 'Japan2012_2020'\n", + "use_auxiliary_plane = False\n", "\n", - " # t indices are 11, 12, 13, 22, 23\n", - " # stress tensor in ENU\n", - " # Tij = np.array([[t[0], t[1], t[2]], [t[1], t[3], t[4]], [t[2], t[4], -t[0]-t[3]]])\n", - " # stress tensor in ESD\n", - " Tij = np.array([[t[0], -t[1], -t[2]], [-t[1], t[3], t[4]], [-t[2], t[4], -t[0] - t[3]]])\n", - " return Tij\n", + "if mydata == \"NZ\":\n", + " df = pd.read_csv(\"GeoNet_CMT_solutions.csv\")\n", + " df = df.rename(columns={\"Latitude\": \"lat\", \"Longitude\": \"lon\", \"CD\": \"depth\"})\n", + " # Kaikoura\n", + " lon0, lon1, lat0, lat1, depth0, depth1 = 172.7, 174.9, -42.7, -40.5, 0, 20\n", + " # region around the Alpine Fault\n", + " # lon0, lon1, lat0, lat1, depth0, depth1 = 166.6, 170.6, -44.6, -43.3, 0, 20\n", + " # North of Northern island\n", + " # lon0, lon1, lat0, lat1, depth0, depth1 = 176.1, 178.9, -40.3, -36.8, 15, np.inf\n", + " # South of Southern island\n", + " # lon0, lon1, lat0, lat1, depth0, depth1 = 163.1, 169.5, -49.1, -44.9, 15, np.inf\n", "\n", + "elif mydata == \"Japan2000_2010\":\n", + " # Downloaded from http://www.isc.ac.uk/iscbulletin/search/fmechanisms/#reviewed\n", + " df = pd.read_csv(mydata + \".csv\", skiprows=23 if mydata == \"Japan2000_2010\" else 25)\n", + " df = df.rename(columns={\"LAT \": \"lat\", \"LON \": \"lon\", \"DEPTH\": \"depth\"})\n", + " df = df.rename(columns={\"STRIKE\": \"strike1\", \"DIP \": \"dip1\", \"RAKE \": \"rake1\"})\n", + " df = df.rename(columns={\"STRIKE.1\": \"strike2\", \"DIP .1\": \"dip2\", \"RAKE .1\": \"rake2\"})\n", + " lon0, lon1, lat0, lat1, depth0, depth1 = 141, np.inf, -np.inf, 38.2, 20, np.inf\n", + "elif mydata == \"\":\n", + " # default file from Stressinverse_1.1.2\n", + " print(\"using default dataset\")\n", + " strike, dip, rake = np.loadtxt(\"West_Bohemia_mechanisms.dat\", skiprows=2, unpack=True)\n", + "else:\n", + " raise ValueError(\"unknown dataset\")\n", + "\n", + "if mydata != \"\":\n", + " df = df[[\"lat\", \"lon\", \"depth\", \"strike1\", \"dip1\", \"rake1\", \"strike2\", \"dip2\", \"rake2\"]]\n", + " df = df.apply(pd.to_numeric, errors=\"coerce\")\n", + " df.dropna(inplace=True)\n", + "\n", + " df = df[\n", + " (df[\"lon\"] > lon0)\n", + " & (df[\"lon\"] < lon1)\n", + " & (df[\"lat\"] > lat0)\n", + " & (df[\"lat\"] < lat1)\n", + " & (df[\"depth\"] > depth0)\n", + " & (df[\"depth\"] < depth1)\n", + " ]\n", + " print(df)\n", + "\n", + " strike, dip, rake = [df[name].values for name in [\"strike1\", \"dip1\", \"rake1\"]]\n", + " strikeAP, dipAP, rakeAP = [df[name].values for name in [\"strike2\", \"dip2\", \"rake2\"]]\n", + "print(strike)\n", "\n", - "def ComputeSortedEigenValuesAndVectors(Tij):\n", - " \"\"\"Compute eigen Values and eigen Vectors, to compute R ratio and SH_max\"\"\"\n", + "if use_auxiliary_plane:\n", + " # set for 30% of focal mechanism the wrong plane\n", + " nfocal = len(strike)\n", + " listn = np.random.choice(range(nfocal), size=int(0.3 * nfocal), replace=False)\n", + " strike[listn] = strikeAP[listn]\n", + " dip[listn] = dipAP[listn]\n", + " rake[listn] = rakeAP[listn]\n", "\n", - " # each column of eigenVectors is an eigenvector\n", - " eigenValues, eigenVectors = np.linalg.eig(Tij)\n", - " # important notes:\n", - " # the largest principal stress corresponds with the smallest seigeValue!\n", - " # the eigenValues should not be sorted in absolute value!\n", - " idx = eigenValues.argsort()\n", - " eigenValues = eigenValues[idx]\n", - " eigenVectors = eigenVectors[:, idx]\n", - " return eigenValues, eigenVectors" + "strike0 = np.radians(strike)\n", + "dip0 = np.radians(dip)\n", + "rake0 = np.radians(rake)" ] }, { @@ -488,78 +495,23 @@ " S = ComputeSlipVector(strike, dip, rake)\n", " A = AssembleAmatrix(FaultNormal)\n", " Tij = computeStressTensor(A, S)\n", - " eigenValues[irealization, :], eigenVectors[irealization, :, :] = ComputeSortedEigenValuesAndVectors(Tij)\n", - "\n", - "# print('stress tensor in ESD\\n', Tij)\n", - "# print(\"sorted eigen values\", eigenValues)\n", - "# print(\"sorted eigen vectors\\n\", eigenVectors)\n", + " (\n", + " eigenValues[irealization, :],\n", + " eigenVectors[irealization, :, :],\n", + " ) = ComputeSortedEigenValuesAndVectors(Tij)\n", "\n", "# plot the eigen vector in the focal sphere\n", "v1 = (eigenVectors[:, :, 0]).T\n", "v2 = (eigenVectors[:, :, 1]).T\n", "v3 = (eigenVectors[:, :, 2]).T\n", - "plot_stress_axes(v1, v2, v3, \"\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def ComputeSHmax(eigenVectors):\n", - " # largest principal stress\n", - " uz = np.array([0.0, 0.0, 1.0])\n", - " idv = int(np.argmax(np.absolute(np.dot(eigenVectors.T, uz))))\n", - " # print(\"most vertical vector\", idv + 1)\n", - "\n", - " # compute largest eigen value not the most vertical eigen vector\n", - " if idv == 0:\n", - " ih = 1\n", - " else:\n", - " ih = 0\n", - "\n", - " vh = eigenVectors[:, ih]\n", - " # SHmax: azimuth of the maximum horizontal compressive stress\n", - " SHmax = np.degrees(np.arctan2(vh[1], vh[0]))\n", - " return SHmax\n", - "\n", - "\n", - "def plotSHmaxR(eigenVectors, eigenValues):\n", - " # Compute SH_max and R (s2_ratio)\n", - " R = (eigenValues[:, 0] - eigenValues[:, 1]) / (eigenValues[:, 0] - eigenValues[:, 2])\n", - " plt.title(\"R distribution\", fontsize=16)\n", - " plt.hist(R)\n", - " plt.show()\n", - " nobservations = eigenValues.shape[0]\n", - " SHmax = np.zeros(nobservations)\n", - " for i in range(nobservations):\n", - " SHmax[i] = ComputeSHmax((eigenVectors[i, :, :]).T)\n", - " plt.title(r\"$SH_{max}$ distribution\", fontsize=16)\n", - " plt.hist(SHmax)\n", - " plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "plot_stress_axes(v1, v2, v3, \"\")\n", "plotSHmaxR(eigenVectors, eigenValues)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -573,7 +525,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.10" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/requirements.txt b/requirements.txt index 4481153..468d805 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,2 @@ -numpy==1.16.* -matplotlib==3.* -seaborn==0.11.* -scipy==1.5.* +numpy==1.21.5 +pandas==1.4.1