forked from jolivetr/csi
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathCDM.py
222 lines (175 loc) · 8.16 KB
/
CDM.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
'''
A Compound Dislocation Model (CDM) sub-class.
Written by T. Shreve, June 2019
Manual cleaned up by R. Jolivet, Feb 2024
'''
# Import Externals stuff
import numpy as np
from argparse import Namespace
import warnings
# Personals
from . import CDMfull
from .Pressure import Pressure
#sub-class CDM
class CDM(Pressure):
# ----------------------------------------------------------------------
# Initialize class
def __init__(self, name, utmzone=None, ellps='WGS84', lon0=None, lat0=None,
verbose=True):
'''
Sub-class implementing CDM pressure objects.
Args:
* name : Name of the pressure source.
* utmzone : UTM zone (optional, default=None)
Kwargs:
* x0, y0 : Center of pressure source in lat/lon or utm
* z0 : Depth
* ax, ay, az : Semi-axes of the CDM along the x, y and z axes respectively, before applying rotations.
* dip : Clockwise around N-S (Y) axis; dip = 90 means vertically elongated source
* strike : Clockwise from N; strike = 0 means source is oriented N-S
* plunge : Clockwise along E-W (X) axis
* ellps : ellipsoid (optional, default='WGS84')
'''
warnings.warn("CDM is not yet fully tested as the amplitudes are awkward. Please check the code.", FutureWarning)
# Attributes only for CDM and pCDM
# Potency given semi-axes and unit opening
self.DV = None
# Final potency change scaled by self.DV
self.deltapotency = None
# Tensile component on each rectangular dislocation
self.deltaopening = None
# Base class init
super(CDM,self).__init__(name, utmzone=utmzone, ellps=ellps, lon0=lon0, lat0=lat0, verbose=verbose)
self.verbose=verbose
self.source = "CDM"
return
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# Find unit potency change (opening=1)
def computePotency(self):
'''
Computes potency change (m3) of three orthogonal point tensile dislocations, given the semimajor axes.
Returns:
* A : Ratio of horizontal potency (volume) to total potency (volume) variation
* B : Ratio of vertical potency (volume) variation
* self.DV : Change in potency
'''
DVx = 4*self.ellipshape['ay']*self.ellipshape['az'] #unit potency (assume opening = 1 m )
DVy = 4*self.ellipshape['ax']*self.ellipshape['az']
DVz = 4*self.ellipshape['ax']*self.ellipshape['ay']
A = DVz / (DVx + DVy + DVz)
B = DVy / (DVx + DVy)
self.DV = (DVx + DVy + DVz)
return A, B, self.DV
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# Some building routines that can be touched... I guess
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# Set volume and pressure changes
def setOpening(self, opening):
'''
Set deltavolume given deltapressure.
Returns:
* deltapressure : Pressure change
'''
self.deltaopening = opening
self.opening2potency()
# All done
return
# ----------------------------------------------------------------------
#
# ----------------------------------------------------------------------
# Convert openign to change in potency
def opening2potency(self):
'''
Converts opening (m) to potency change (m3) for CDM.
Uses formula from :
Nikkhoo, M., Walter, T. R., Lundgren, P. R., and Prats-Iraola, P. (2017). Compound dislocation models(CDMs) for volcano deformation analyses.Geophysical Journal International, 208(2):877–894
Returns:
* deltapotency : strength of dislocation source, or volume available to host fluids intruding into cavity from the outside
'''
ax, ay, az = self.ellipshape['ax'], self.ellipshape['ay'], self.ellipshape['az']
self.deltapotency = (4*(ax*ay + ay*az + ax*az)*self.deltaopening)
# All done
return self.deltapotency
def pressure2dis(self, data, delta="volume", volume=None):
'''
Computes the surface displacement at the data location using pCDM. ~~~ This is where the good stuff happens ~~
Args:
* data : data object from gps or insar.
* delta : total potency change, in units of volume
Returns:
* u : x, y, and z displacements
'''
# Get parameters
if self.ellipshape is None:
raise Exception("Error: Need to define shape of source (run self.createShape)")
#define dictionary entries as variables
ellipse = Namespace(**self.ellipshape)
# Set the potency change
if volume is None:
self.DV = (4*(ellipse.ax*ellipse.ay + ellipse.ax*ellipse.az + ellipse.ay*ellipse.az))
opening = self.deltapotency/self.DV
print("Opening is", opening)
else:
# Set the potency value
if delta=="volume":
opening = 1.0
self.DV = (4*(ellipse.ax*ellipse.ay + ellipse.ax*ellipse.az + ellipse.ay*ellipse.az))
if delta=="pressure":
raise Exception("Error: CDM requires potencies, not pressure")
# Set poisson's ratio
if self.nu is None:
self.nu = 0.25
# Get data position -- in m
x = data.x*1000
y = data.y*1000
z = np.zeros(x.shape) # Data are at the surface
# Run it for pCDM pressure source
# Convert from dip, strike, and plunge to omegaX, omegaY, omegaZ - verify they have same convention as Yang
omegaX = ellipse.plunge
omegaY = 90. - ellipse.dip
omegaZ = ellipse.strike
if (opening!=0.0):
#### x0, y0 need to be in utm and meters
Ux,Uy,Uz,Dv = CDMfull.displacement(x, y, z, ellipse.x0m*1000, ellipse.y0m*1000, ellipse.z0, omegaX, omegaY, omegaZ, ellipse.ax, ellipse.ay, ellipse.az, opening, self.nu)
else:
dp_dis = np.zeros((len(x), 3))
# All done
# concatenate into one array
u = np.vstack([Ux, Uy, Uz]).T
return u
# ----------------------------------------------------------------------
# Define the shape of the CDM
#
# ----------------------------------------------------------------------
def createShape(self, x, y, z0, ax, ay, az, dip, strike, plunge, latlon=True):
""""
Defines the shape of the CDM pressure source.
Args:
* x0, y0 : Center of pressure source in lat/lon or utm
* z0 : Depth
* ax, ay, az : Semi-axes of the CDM along the x, y and z axes respectively, before applying rotations.
* dip : Clockwise around N-S (Y) axis; dip = 90 means vertical source
* strike : Clockwise from N; strike = 0 means source is oriented N-S
* plunge : Clockwise along E-W (X) axis
Returns:
None
"""
if latlon is True:
self.lon, self.lat = x, y
lon1, lat1 = x, y
self.pressure2xy()
x1, y1 = self.xf, self.yf
else:
self.xf, self.yf = x, y
x1, y1 = x, y
self.pressure2ll()
lon1, lat1 = self.lon, self.lat
self.ellipshape = {'x0': lon1,'x0m': x1,'y0': lat1, 'y0m': y1, 'z0': z0,'ax': ax,'ay': ay, 'az': az, 'dip': dip,'strike': strike, 'plunge': plunge}
return
#EOF