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)