Skip to content

Commit

Permalink
added dispersion for 1D models and colors
Browse files Browse the repository at this point in the history
  • Loading branch information
@pmoulik
pmoulik committed Aug 13, 2024
1 parent 0110866 commit c014fdc
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 27 deletions.
10 changes: 10 additions & 0 deletions avni/constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,16 @@
ureg = None # this is initialized to MKS units in _init_.py
customunits = 'units.ini'

# custome colors for color-blindness
# use these instead of RGB, also prints well in B/Warning
colors = {
'black':'#231F20',
'yellow':'#FFBB00',
'red':'#FB6542',
'blue':'#375E97'
}


#Color scales
colorscale = {
'avni': {'name': 'avni','description': \
Expand Down
93 changes: 86 additions & 7 deletions avni/models/reference1d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
import warnings
import os
import pdb
import typing as tp

if sys.version_info[0] >= 3: unicode = str

Expand Down Expand Up @@ -564,8 +565,9 @@ def get_Love_elastic(self):

# Add metadata
for field in ['a','c','n','l','f','vp','vs','vphi','xi','phi','Zp', 'Zs','kappa','mu','as','ap']:
self.metadata['parameters'].append(field)
self.metadata['units'].append(str(self.data[field].pint.units))
if field not in self.metadata['parameters']:
self.metadata['parameters'].append(field)
self.metadata['units'].append(str(self.data[field].pint.units))

def get_mineralogical(self):
'''
Expand Down Expand Up @@ -612,14 +614,16 @@ def get_mineralogical(self):

# Add metadata
for field in ['gravity','Brunt-Vaisala','Bullen','pressure']:
self.metadata['parameters'].append(field)
self.metadata['units'].append(str(self.data[field].pint.units))
if field not in self.metadata['parameters']:
self.metadata['parameters'].append(field)
self.metadata['units'].append(str(self.data[field].pint.units))

# Poisson ratio
vsbyvp1d = self.data['vs'].div(self.data['vp'])
self.data['poisson2'] = np.divide((1.-(2.*vsbyvp1d*vsbyvp1d)),(2.-(2.*vsbyvp1d*vsbyvp1d)))
self.metadata['parameters'].append('poisson')
self.metadata['units'].append('dimensionless')
self.data['poisson'] = np.divide((1.-(2.*vsbyvp1d*vsbyvp1d)),(2.-(2.*vsbyvp1d*vsbyvp1d)))
if 'poisson' not in self.metadata['parameters']:
self.metadata['parameters'].append('poisson')
self.metadata['units'].append('dimensionless')

# delete a file
with contextlib.suppress(FileNotFoundError): os.remove(file)
Expand Down Expand Up @@ -718,6 +722,81 @@ def get_discontinuity(self):
if itopmantle >0: disc['itopmantle'] = itopmantle-1
self.metadata['discontinuities'] = disc

def get_dispersion(self,
period: tp.Union[int,float,np.int64]):
'''
Computes the dispersion correction factors and updates the elastic parameters
to those valid at the requested period.
'''
import pint_pandas

ref_period = self.metadata['ref_period']
if ref_period != 1: raise NotImplementedError('Dispersion cannot be implemented when reference period does not equal 1 second for the reported parameters.')

if self.data is None or self._nlayers == 0:
warnings.warn('reference1D data arrays are not allocated. Trying to apply coefficients_to_cards from within get_dispersion')
self.coefficients_to_cards()

for param in ['a','c','n','l','f']:
if param not in self.data.keys():
warnings.warn('Love parameter '+param+' is not allocated. Trying to apply get_Love_elastic from within get_dispersion')
self.get_Love_elastic()

A = self.data['a'].pint.to_base_units().values.quantity.magnitude
C = self.data['c'].pint.to_base_units().values.quantity.magnitude
N = self.data['n'].pint.to_base_units().values.quantity.magnitude
L = self.data['l'].pint.to_base_units().values.quantity.magnitude
F = self.data['f'].pint.to_base_units().values.quantity.magnitude
rho = self.data['rho'].pint.to_base_units().values.quantity.magnitude
qka = ((1/self.data['qkappa']).replace(np.inf, 0)).values.quantity.magnitude
qmu = ((1/self.data['qmu']).replace(np.inf, 0)).values.quantity.magnitude

# correction factors
E = (4/3*(A+C-2*F+5*N+6*L)/(8*A+3*C+4*F+8*L))
xac = 1 - 2*np.log(period)/np.pi*((1-E)*qka+E*qmu)
xnl = 1 - 2*np.log(period)/np.pi*qmu
xf = 1 - 2*np.log(period)/np.pi*(((1-E)*qka-E/2*qmu)/(1-3/2*E))
xvp = 1 - np.log(period)/np.pi*((1-E)*qka+E*qmu)
xvs = 1 - np.log(period)/np.pi*qmu

# Apply the correction
A *= xac
C *= xac
N *= xnl
L *= xnl
F *= xf

# Updte data fields
PA_ = pint_pandas.PintArray; temp_dict = {}
temp_dict['vph'] = PA_(np.sqrt(A/rho),dtype="pint[m/s]")
temp_dict['vpv'] = PA_(np.sqrt(C/rho),dtype="pint[m/s]")
temp_dict['vsh'] = PA_(np.sqrt(N/rho),dtype="pint[m/s]")
temp_dict['vsv'] = PA_(np.sqrt(L/rho),dtype="pint[m/s]")
temp_dict['eta'] = PA_(F/(A-2*L),dtype="pint[dimensionless]")
temp_dict['E'] = PA_(E,dtype="pint[dimensionless]")
temp_dict['xac'] = PA_(xac,dtype="pint[dimensionless]")
temp_dict['xnl'] = PA_(xnl,dtype="pint[dimensionless]")
temp_dict['xf'] = PA_(xf,dtype="pint[dimensionless]")
temp_dict['xvp'] = PA_(xvp,dtype="pint[dimensionless]")
temp_dict['xvs'] = PA_(xvs,dtype="pint[dimensionless]")

modelarr = pd.DataFrame(temp_dict)
self.data.drop(columns=['vph','vpv','vsh','vsv','eta'],inplace=True)
for field in ['E','xac','xnl','xf','xvp','xvs']:
if field in self.data.columns: self.data.drop(columns=[field],inplace=True)
self.data = pd.concat([self.data,modelarr],axis=1)
self.metadata['ref_period'] = period

# redo calculations
self.get_Love_elastic()
self.get_discontinuity()
getmineral = False
for field in ['gravity','Brunt-Vaisala','Bullen','pressure','poisson']:
if field in self.data.columns:
self.data.drop(columns=[field])
getmineral = True
if getmineral: self.get_mineralogical()

def get_custom_parameter(self,parameters):
'''
Get the arrays of custom parameters defined in various Earth models
Expand Down
41 changes: 21 additions & 20 deletions avni/plots/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
from .. import tools
from .. import data
from .. import constants
from ..constants import colors
from .common import initializecolor,updatefont

##########################################################################
Expand Down Expand Up @@ -1498,13 +1499,13 @@ def plotreference1d(ref1d,
gs = gridspec.GridSpec(3, 1, height_ratios=height_ratios)
fig.patch.set_facecolor('white')
ax01=plt.subplot(gs[0])
ax01.plot(depthkmarr,rho,'k')
ax01.plot(depthkmarr,vsv,'b')
ax01.plot(depthkmarr,vsh,'b:')
ax01.plot(depthkmarr,vpv,'r')
ax01.plot(depthkmarr,vph,'r:')
ax01.plot(depthkmarr,rho,'-',color=colors['black'])
ax01.plot(depthkmarr,vsv,'-',color=colors['blue'])
ax01.plot(depthkmarr,vsh,':',color=colors['blue'])
ax01.plot(depthkmarr,vpv,'-',color=colors['red'])
ax01.plot(depthkmarr,vph,':',color=colors['red'])
mantle=np.where( depthkmarr < 2891.)
ax01.plot(depthkmarr[mantle],eta[mantle],'g')
ax01.plot(depthkmarr[mantle],eta[mantle],colors['yellow'])
ax01.set_xlim([0., constants.R.to('km').magnitude])
ax01.set_ylim([0, 14])

Expand All @@ -1525,7 +1526,7 @@ def plotreference1d(ref1d,
ax01.xaxis.set_minor_locator(minorLocator)
ax01.set_ylabel('Velocity (km/sec), density (g/cm'+'$^3$'+') or '+'$\eta$')

for para,color,xloc,yloc in [("$\eta$",'g',1500.,2.),("$V_S$",'b',1500.,7.8),("$V_P$",'r',1500.,13.5),("$\\rho$",'k',1500.,4.5),("$V_P$",'r',4000.,9.2),("$\\rho$",'k',4000.,12.5),("$V_S$",'b',5500.,4.5)]:
for para,color,xloc,yloc in [("$\eta$",colors['yellow'],1500.,2.),("$V_S$",colors['blue'],1500.,7.8),("$V_P$",colors['red'],1500.,13.5),("$\\rho$",colors['black'],1500.,4.5),("$V_P$",colors['red'],4000.,9.2),("$\\rho$",colors['black'],4000.,12.5),("$V_S$",colors['blue'],5500.,4.5)]:
ax01.annotate(para,color=color,
xy=(3, 1), xycoords='data',
xytext=(xloc/constants.R.to('km').magnitude, yloc/14.), textcoords='axes fraction',
Expand All @@ -1534,27 +1535,27 @@ def plotreference1d(ref1d,

ax11=plt.subplot(gs[1])
depthselect=np.intersect1d(np.where( depthkmarr >= zoomdepth[0]),np.where( depthkmarr <= zoomdepth[1]))
ax11.plot(depthkmarr[depthselect],rho[depthselect],'k')
ax11.plot(depthkmarr[depthselect],rho[depthselect],color=colors['black'])
if isotropic:
ax11.plot(depthkmarr[depthselect],vs[depthselect],'b')
ax11.plot(depthkmarr[depthselect],vs[depthselect],color=colors['blue'])
else:
ax11.plot(depthkmarr[depthselect],vsv[depthselect],'b')
ax11.plot(depthkmarr[depthselect],vsh[depthselect],'b:')
ax11.plot(depthkmarr[depthselect],vsv[depthselect],color=colors['blue'])
ax11.plot(depthkmarr[depthselect],vsh[depthselect],':',color=colors['blue'])
ax12 = ax11.twinx()
if isotropic:
ax12.plot(depthkmarr[depthselect],vp[depthselect],'r')
ax12.plot(depthkmarr[depthselect],vp[depthselect],color=colors['red'])
else:
ax12.plot(depthkmarr[depthselect],vpv[depthselect],'r')
ax12.plot(depthkmarr[depthselect],vph[depthselect],'r:')
ax12.plot(depthkmarr[depthselect],vpv[depthselect],color=colors['red'])
ax12.plot(depthkmarr[depthselect],vph[depthselect],':',color=colors['red'])

ax11.plot(depthkmarr[depthselect],eta[depthselect],'g')
ax11.plot(depthkmarr[depthselect],eta[depthselect],color=colors['yellow'])
ax11.set_xlim(zoomdepth)
ax11.set_ylim([0, 7])
ax12.set_xlim(zoomdepth)
ax12.set_ylim([-2, 12])
ax11.set_ylabel('Shear velocity (km/sec), density (g/cm'+'$^3$'+') or '+'$\eta$')
ax12.set_ylabel('Compressional velocity (km/sec)')
for para,color,xloc,yloc in [("$\eta$",'g',150.,1.),("$V_{S}$",'b',150.,4.3),("$V_{P}$",'r',120.,5.5),("$\\rho$",'k',150.,3.8)]:
for para,color,xloc,yloc in [("$\eta$",colors['yellow'],150.,1.),("$V_{S}$",colors['blue'],150.,4.3),("$V_{P}$",colors['red'],120.,5.5),("$\\rho$",colors['black'],150.,3.8)]:
ax11.annotate(para,color=color,
xy=(3, 1), xycoords='data',
xytext=(xloc/1000., yloc/7.), textcoords='axes fraction',
Expand All @@ -1570,8 +1571,8 @@ def plotreference1d(ref1d,


ax21=plt.subplot(gs[2], sharex=ax11)
ax21.plot(depthkmarr[depthselect],anisoVs[depthselect],'b')
ax21.plot(depthkmarr[depthselect],anisoVp[depthselect],'r')
ax21.plot(depthkmarr[depthselect],anisoVs[depthselect],color=colors['blue'])
ax21.plot(depthkmarr[depthselect],anisoVp[depthselect],color=colors['red'])
ax21.set_ylim([0, 5])
ax21.set_xlim(zoomdepth)
majorLocator = MultipleLocator(1)
Expand All @@ -1581,15 +1582,15 @@ def plotreference1d(ref1d,
ax21.yaxis.set_major_formatter(majorFormatter)
# for the minor ticks, use no labels; default NullFormatter
ax21.yaxis.set_minor_locator(minorLocator)
for para,color,xloc,yloc in [('Q'+'$_{\mu}$','k',400.,2.5),("$a_{S}$",'b',150.,3.7),("$a_{P}$",'r',50.,3.5)]:
for para,color,xloc,yloc in [('Q'+'$_{\mu}$',colors['black'],400.,2.5),("$a_{S}$",colors['blue'],150.,3.7),("$a_{P}$",colors['red'],50.,3.5)]:
ax21.annotate(para,color=color,
xy=(3, 1), xycoords='data',
xytext=(xloc/1000., yloc/4.), textcoords='axes fraction',
horizontalalignment='left', verticalalignment='top')


ax22 = ax21.twinx()
ax22.plot(depthkmarr[depthselect],qmu[depthselect],'k')
ax22.plot(depthkmarr[depthselect],qmu[depthselect],color=colors['black'])
ax21.set_xlabel('Depth (km)')
ax21.set_ylabel("$V_P$"+' or '+"$V_S$"+' anisotropy (%)')
ax22.set_ylabel('Shear attenuation Q'+'$_{\mu}$')
Expand Down

0 comments on commit c014fdc

Please sign in to comment.