diff --git a/avni/constants.py b/avni/constants.py index e888e41..8222998 100644 --- a/avni/constants.py +++ b/avni/constants.py @@ -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': \ diff --git a/avni/models/reference1d.py b/avni/models/reference1d.py index c4b1d5f..8e79b1b 100644 --- a/avni/models/reference1d.py +++ b/avni/models/reference1d.py @@ -26,6 +26,7 @@ import warnings import os import pdb +import typing as tp if sys.version_info[0] >= 3: unicode = str @@ -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): ''' @@ -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) @@ -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 diff --git a/avni/plots/models.py b/avni/plots/models.py index 9c8ac52..5f6c9ee 100644 --- a/avni/plots/models.py +++ b/avni/plots/models.py @@ -45,6 +45,7 @@ from .. import tools from .. import data from .. import constants +from ..constants import colors from .common import initializecolor,updatefont ########################################################################## @@ -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]) @@ -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', @@ -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', @@ -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) @@ -1581,7 +1582,7 @@ 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', @@ -1589,7 +1590,7 @@ def plotreference1d(ref1d, 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}$')