subfaults.py.html
Source file:   subfaults.py
Directory:   /var/www/html/clawpack/links/awr11/chile2010awr
Converted:   Sat Feb 5 2011 at 15:39:51   using clawcode2html
This documentation file will not reflect any later changes in the source file.

 

from pylab import *
from pyclaw.geotools.okada import okadamap
import re

Rearth = 6367.5e3
class FaultModel(object):
    pass

def read_subfault_model(fname):
    lines = open(fname,'r').readlines()
    regexp_dx = re.compile(r"Dx=[ ]*(?P[^k]*)")
    regexp_dy = re.compile(r"Dy=[ ]*(?P[^k]*)")
    regexp_nx = re.compile(r"nx[^=]*=[ ]*(?P[^D]*)")
    regexp_ny = re.compile(r"ny[^=]*=[ ]*(?P[^D]*)")
    
    i_dxdy = -1
    for i,line in enumerate(lines):
        result_dx = regexp_dx.search(line)
        result_dy = regexp_dy.search(line)
        result_nx = regexp_nx.search(line)
        result_ny = regexp_ny.search(line)
        if result_dx and result_dy:
            i_dxdy = i
            dx = float(result_dx.group('dx')) * 1.e3
            dy = float(result_dy.group('dy')) * 1.e3
            nx = int(result_nx.group('nx'))
            ny = int(result_ny.group('ny'))
            print "Subfaults: dx = %s, dy = %s, nx = %s, ny = %s " \
                  % (dx, dy, nx, ny)
            
    if i_dxdy == -1:
        print "*** Did not find a line containing both Dx and Dy"
        raise Exception()
    
    fault_bdry = []
    for i in range(i_dxdy+3, i_dxdy+8):    
        fault_bdry.append([float(s) for s in lines[i].split()])
    fault_bdry = array(fault_bdry)
    #print "+++ fault_bdry = ", fault_bdry

        
    
    data = loadtxt(fname,skiprows=9+i_dxdy)
    arrayshape = (ny,nx)
    fm = FaultModel()
    fm.arrayshape = arrayshape
    fm.subfault_length = dx
    fm.subfault_width = dy
    fm.hypo_longitude = 0.25*(sum(fault_bdry[1:,0]))
    fm.hypo_latitude = 0.25*(sum(fault_bdry[1:,1]))
    
    fm.latitude = reshape(data[:,0],arrayshape)
    fm.longitude = reshape(data[:,1],arrayshape)
    fm.depth = reshape(data[:,2],arrayshape) * 1.e3 # convert km to meters
    fm.slip = reshape(data[:,3],arrayshape) * 0.01  # convert cm to meters
    fm.rake = reshape(data[:,4],arrayshape)
    fm.strike = reshape(data[:,5],arrayshape)
    fm.dip = reshape(data[:,6],arrayshape)
    fm.fault_bdry = fault_bdry
    
    return fm
    
def make_okada_dz(fm, faultparams):
    mx = faultparams['mx']
    my = faultparams['my']
    
    X=linspace(faultparams['xlower'],faultparams['xupper'],mx)
    Y=linspace(faultparams['ylower'],faultparams['yupper'],my)
    dZ = zeros((my,mx))
    
    okadaparams = {}
    
    for i in range(fm.arrayshape[0]):
        for j in range(fm.arrayshape[1]): 
            okadaparams["Focal_Depth"] = fm.depth[i,j]
            okadaparams["Fault_Length"] = fm.subfault_length
            okadaparams["Fault_Width"] = fm.subfault_width
            okadaparams["Dislocation"] = fm.slip[i,j]
            okadaparams["Strike_Direction"] = fm.strike[i,j]
            okadaparams["Dip_Angle"] = fm.dip[i,j]
            okadaparams["Slip_Angle"] = fm.rake[i,j]
            okadaparams["Epicenter_Latitude"] = fm.latitude[i,j]
            okadaparams["Epicenter_Longitude"] = fm.longitude[i,j]
            
            dZ = dZ + okadamap(okadaparams, X, Y)
            
    return X,Y,dZ
    
    
def write_dz(fname,X,Y,dZ):
    fid = open(fname, 'w')
    t0 = 0.
    tend = 100.
    nt = 2
    T = linspace(t0, tend, nt)
    for it in T:
        alpha=(it-t0)/(tend-t0)
        
        for jj in range(len(Y)):
            j=-1-jj
            for i in range(len(X)) :
                fid.write('%012.6e %012.6e %012.6e %012.6e \n' % (it,X[i],Y[j],alpha*dZ[j,i]))

    fid.close()
    print "Created ",fname
    
    
def make_dz_chile2010():
    fname_subfaults = 'usgs_subfault.txt'
    fname_dtopo = 'usgs_subfault.tt1'
    faultparams = {}
    faultparams['mx'] = 120
    faultparams['my'] = 160
    faultparams['xlower'] = -76
    faultparams['xupper'] = -70
    faultparams['ylower'] = -40
    faultparams['yupper'] = -32
    fm = read_subfault_model(fname_subfaults)
    X,Y,dZ = make_okada_dz(fm, faultparams)
    write_dz(fname_dtopo, X,Y,dZ)
    return X,Y,dZ


def plot_dz(X,Y,dZ):
    from pyclaw.plotters import colormaps
    figure(200)
    clf()
    bwr = colormaps.blue_white_red
    pcolor(X,Y,dZ,cmap=bwr)
    clim(-3,3)
    axis('scaled')
    colorbar()
    
if __name__=="__main__":
    X,Y,dZ = make_dz_chile2010()
    #plot_dz(X,Y,dZ)