In [1]:
import linecache
import numpy as np
In [2]:
print("# numpy version = ",np.__version__)

# numpy version =  1.16.4
In [3]:
def rho_from_vp_ludwig(vp):
    """
    Brocher (2005) equation (1)
    """
    rho = 1.6612*vp
    rho += -0.4721*(vp**2)
    rho += 0.0671*(vp**3)
    rho += -0.0043*(vp**4)
    rho += 0.000106*(vp**5)
    return rho

def rho_from_vp_gardner(vp):
    rho = 1.74*(vp**0.25)
    return rho

def vp_from_rho(rho):
    """
    Brocher (2005) equation (5)
    """
    vp = 39.128*rho
    vp += -63.064*(rho**2)
    vp += 37.083*(rho**3)
    vp += -9.1819*(rho**4)
    vp += 0.8228*(rho**5)
    return vp

def vs_from_vp(vp):
    """
    Brocher (2005) equation (6)
    """
    vs = 0.7858 
    vs += -1.2344*vp
    vs += 0.7949*(vp**2)
    vs += -0.1238*(vp**3)
    vs += 0.0064*(vp**4)
    return vs


def vp_from_vs(vs):
    """
    Brocher (2005) equation (9)
    """
    vp = 0.9409
    vp += 2.0947*vs
    vp += -0.8206*(vs**2)
    vp += 0.2683*(vs**3)
    vp += -0.0251*(vs**4)
    return vp

def qs_from_vs(vs):
    """
    Brocher (2008) equation (19)
    """
    qs = -16 + 104.13*vs
    qs += -25.225*(vs**2)
    qs += 8.2184*(vs**3)
    return qs
def qp_from_qs(qs):
    return 2*qs

# Testing
#vs = 1.5
#vp = (3**0.5)*vs
#vp = vp_from_vs(vs)
##vp = 0.5
#rho_l = rho_from_vp_ludwig(vp)
#rho_g = rho_from_vp_gardner(vp)
##vs = vs_from_vp(vp)
#print ('vp: ', vp, ' vs: ', vs, ' rho: ', rho_l, rho_g)
#
#print('vp from rho: ', vp_from_rho(rho_l))

#qs = qs_from_vs(vs)
#qp = qp_from_qs(qs)
#print('qs: ', qs, ' qp: ', qp)
In [4]:
# -2 would be needed from Cliff's email. the first and last grid point would have fixed value? i.e., not inveted by data.
# total line number of his file is 858. this support 66*13
# also each line has 79 values not 81.

nx = 81-2 # 79
ny = 68-2 # 66
nz = 15-2 # 13

vs = np.zeros((nx,ny,nz))
vp = np.zeros_like(vs)

f_vs = open('Vs_model.dat', 'r')
f_vp = open('Vp_model.dat', 'r')
In [5]:
vs_file = "Vs_model.dat"
vp_file = "Vp_model.dat"

nx = 81-2 # 79
ny = 68-2 # 66
nz = 15-2 # 13

vs = np.zeros((nx,ny,nz))
vp = np.zeros_like(vs)

n=0

for z in range(nz):
    for y in range(ny):
        n=n+1
        vs_line = (linecache.getline(vs_file ,n)).split()
        vs_num = len(vs_line)
        vp_line = (linecache.getline(vp_file ,n)).split()
        vp_num = len(vp_line)
        
        if vp_num == vs_num:
            #print("do something")
            #for x in range(vs_line):
            for x, vsOUT in enumerate(vs_line):
                #print('x, vs', x, vs)
                vp[x,y,z]=vp_line[x]
                vs[x,y,z]=vs_line[x]                
                #print('# x = ',x,' y =',y, 'z = ',z,' vp = ',vp[x,y,z], 'vs = ',vs[x,y,z])
                
        else:
            print("something wrong!")
            print('y = ',y, 'z = ',z)
In [6]:
xgrid = [235.9000,235.9500,236.0000,236.0500,236.1000,236.1500,236.2000,236.2500,236.3000,236.3500,236.4000,236.4500,236.5000,236.5500,236.6000,236.6500,236.7000,236.7500,236.8000,236.8500,236.9000,236.9500,237.0000,237.0500,237.1000,237.1500,237.2000,237.2500,237.3000,237.3500,237.4000,237.4500,237.5000,237.5500,237.6000,237.6500,237.7000,237.7500,237.8000,237.8500,237.9000,237.9500,238.0000,238.0500,238.1000,238.1500,238.2000,238.2500,238.3000,238.3500,238.4000,238.4500,238.5000,238.5500,238.6000,238.6500,238.7000,238.7500,238.8000,238.8500,238.9000,238.9500,239.0000,239.0500,239.1000,239.1500,239.2000,239.2500,239.3000,239.3500,239.4000,239.4500,239.5000,239.5500,239.6000,239.6500,239.7000,239.7500,239.8000]
len(xgrid)
Out[6]:
79
In [7]:
print(xgrid)

[235.9, 235.95, 236.0, 236.05, 236.1, 236.15, 236.2, 236.25, 236.3, 236.35, 236.4, 236.45, 236.5, 236.55, 236.6, 236.65, 236.7, 236.75, 236.8, 236.85, 236.9, 236.95, 237.0, 237.05, 237.1, 237.15, 237.2, 237.25, 237.3, 237.35, 237.4, 237.45, 237.5, 237.55, 237.6, 237.65, 237.7, 237.75, 237.8, 237.85, 237.9, 237.95, 238.0, 238.05, 238.1, 238.15, 238.2, 238.25, 238.3, 238.35, 238.4, 238.45, 238.5, 238.55, 238.6, 238.65, 238.7, 238.75, 238.8, 238.85, 238.9, 238.95, 239.0, 239.05, 239.1, 239.15, 239.2, 239.25, 239.3, 239.35, 239.4, 239.45, 239.5, 239.55, 239.6, 239.65, 239.7, 239.75, 239.8]
In [8]:
ygrid=[36.1000,36.1500,36.2000,36.2500,36.3000,36.3500,36.4000,36.4500,36.5000,36.5500,36.6000,36.6500,36.7000,36.7500,36.8000,36.8500,36.9000,36.9500,37.0000,37.0500,37.1000,37.1500,37.2000,37.2500,37.3000,37.3500,37.4000,37.4500,37.5000,37.5500,37.6000,37.6500,37.7000,37.7500,37.8000,37.8500,37.9000,37.9500,38.0000,38.0500,38.1000,38.1500,38.2000,38.2500,38.3000,38.3500,38.4000,38.4500,38.5000,38.5500,38.6000,38.6500,38.7000,38.7500,38.8000,38.8500,38.9000,38.9500,39.0000,39.0500,39.1000,39.1500,39.2000,39.2500,39.3000,39.3500]
len(ygrid)
Out[8]:
66
In [9]:
print(ygrid)

[36.1, 36.15, 36.2, 36.25, 36.3, 36.35, 36.4, 36.45, 36.5, 36.55, 36.6, 36.65, 36.7, 36.75, 36.8, 36.85, 36.9, 36.95, 37.0, 37.05, 37.1, 37.15, 37.2, 37.25, 37.3, 37.35, 37.4, 37.45, 37.5, 37.55, 37.6, 37.65, 37.7, 37.75, 37.8, 37.85, 37.9, 37.95, 38.0, 38.05, 38.1, 38.15, 38.2, 38.25, 38.3, 38.35, 38.4, 38.45, 38.5, 38.55, 38.6, 38.65, 38.7, 38.75, 38.8, 38.85, 38.9, 38.95, 39.0, 39.05, 39.1, 39.15, 39.2, 39.25, 39.3, 39.35]
In [10]:
zgrid = [-1.00,0.00,1.00,2.00,4.00,6.00,8.00,10.00,12.00,15.00,20.00,25.00,35.00]
len(zgrid)
Out[10]:
13
In [11]:
print(zgrid)

[-1.0, 0.0, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 35.0]
In [12]:
# sw4/examples/pfile/caucasus_mod3.ppmod
#Caucasus
#0.25
#7 38.00 39.50
#19 44.50 49.00
#30 0.00 161.00
#-99 -99 -99 -99
#.TRUE.

outfi = 'uw_model.ppmod'
dgrid = 0.05 # 0.05 degress. matched to xgrid and ygrid
with open(outfi, 'w') as f:
    f.write('UWBayArea' + "\n")
    f.write(str(dgrid)+ "\n")
    f.write('{:d} {:.2f} {:.2f}'.format(len(ygrid),ygrid[0],ygrid[-1])+ "\n")
    f.write('{:d} {:.2f} {:.2f}'.format(len(xgrid),xgrid[0],xgrid[-1])+ "\n")
    f.write('{:d} {:.2f} {:.2f}'.format(len(zgrid),zgrid[0],zgrid[-1])+ "\n")
    f.write("-99 -99 -99 -99"+ "\n")
    f.write(".TRUE."+ "\n")
f.close()


#38.00 44.50 30
 #1  0.00  4.13  2.23  2.62   1000    700 
 #2  3.00  4.13  2.23  2.62   1000    700 
 #3  5.00  4.81  2.60  2.71   1000    700 
# 4  7.00  5.36  2.90  2.78   1000    700 

with open(outfi, 'a') as f:
    for y, lat in enumerate(ygrid):
        for x, lon in enumerate(xgrid):
            #print('{:.2f} {:.2f} {:d}'.format(lat,lon,len(zgrid)))
            f.write('{:.2f} {:.2f} {:d}'.format(lat,lon,len(zgrid))+ "\n")
            for z, depth in enumerate(zgrid):
                #print(lat,lon,depth)
                #rho_l, rho_g. use rho_l 
                rho_l = rho_from_vp_ludwig(vp[x,y,z])
                qs = qs_from_vs(vs[x,y,z])
                qp = qp_from_qs(qs)
                #print('{:.2f} {:.2f} {:.2f}'.format(vp[x,y,z],vs[x,y,z],rho_l))
                #print(' {:d} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}'.format((z+1),depth,vp[x,y,z],vs[x,y,z],rho_l,qp,qs))
                f.write(' {:d} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}'.format((z+1),depth,vp[x,y,z],vs[x,y,z],rho_l,qp,qs)+ "\n")

                
                
                
f.close()
In [ ]:
In [ ]: