deleted obsolete files
This commit is contained in:
Marcel Paffrath
parent
2dd36379a8
commit
21ffbcabc8
@ -1,107 +0,0 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
import sys
|
||||
from obspy import read
|
||||
from obspy import Stream
|
||||
from obspy import Trace
|
||||
from datetime import datetime
|
||||
import numpy as np
|
||||
|
||||
from pylot.core.active import surveyUtils
|
||||
from pylot.core.active import seismicshot
|
||||
import activeSeismoPick
|
||||
reload(seismicshot)
|
||||
reload(surveyUtils)
|
||||
|
||||
#####################################################################################
|
||||
# parameter definitions:#############################################################
|
||||
#traceslist = list(np.arange(1, 49, 1)) # traces (1-48)
|
||||
#shotlist = list(np.arange(302, 352, 1)) # shot-files(.dat) (Paffrath: 302-352) (Hauburg: 353-401) (arange+1!)
|
||||
cutwindow = (0, 0.2) # cut out a part of the trace [seconds]
|
||||
tmovwind = 0.3 # size of the moving window
|
||||
windowsize = (5, 0) # windowsize for AIC picker (indices around HOS picks!!)
|
||||
pickwindow = cutwindow # for local max and threshold picker: fraction of the seismogram used (0...1) TO BE DONE: depends on cutwindow!!!!
|
||||
folm = 0.6
|
||||
|
||||
rockeskyll = False
|
||||
if rockeskyll == True:
|
||||
receiverfile = "Geophone_interpoliert_rockes"
|
||||
sourcefile = "Schusspunkte_rockes"
|
||||
obsdir = "/rscratch/minos22/marcel/flachseismik/rockeskyll_200615_270615/"
|
||||
filename = 'survey_rockes.pickle'
|
||||
else:
|
||||
receiverfile = "Geophone_interpoliert_GZB"
|
||||
sourcefile = "Schusspunkte_GZB"
|
||||
obsdir = "/rscratch/minos22/marcel/flachseismik/GZB_26_06_15_01/"
|
||||
filename = 'survey_GZB.pickle'
|
||||
|
||||
# SNR
|
||||
tsignal = 0.03
|
||||
tgap = 0.0007
|
||||
snrthreshold = 1
|
||||
######################################################################################
|
||||
|
||||
vmin = 333
|
||||
vmax = 5500
|
||||
distBinSize = 2
|
||||
|
||||
###########################################
|
||||
############# Settings: ###################
|
||||
thresholdpick=True
|
||||
localmaxpick=False
|
||||
|
||||
if thresholdpick == True: pickmethod = "threshold"
|
||||
if localmaxpick == True: pickmethod = "localmax"
|
||||
|
||||
HosAic = 'hos' # can be 'hos' or 'aic'
|
||||
###########################################
|
||||
|
||||
starttime = datetime.now()
|
||||
|
||||
print '\n--------------------Starting Script at %s -------------------\n' %starttime.time()
|
||||
if thresholdpick == True: print 'Using treshold pickmethod!\n'
|
||||
elif localmaxpick == True: print 'Using local maximum pick method!\n'
|
||||
if HosAic == 'aic': print 'picking with AIC\n'
|
||||
if HosAic == 'hos': print 'picking with HOS\n'
|
||||
|
||||
survey = activeSeismoPick.Survey(obsdir, sourcefile, receiverfile, True)
|
||||
surveyUtils.setFittedSNR(survey.getShotDict())
|
||||
print '\nDone after %s seconds!\n------------------------------------------------------------------------------\n' % (datetime.now() - starttime).seconds
|
||||
|
||||
count = 0; tpicksum = starttime - starttime
|
||||
|
||||
for shot in survey.data.values():
|
||||
tstartpick = datetime.now(); count += 1
|
||||
for traceID in shot.getTraceIDlist():
|
||||
distance = shot.getDistance(traceID) # receive distance
|
||||
|
||||
pickwin_used = pickwindow # use pickwindow set in the parameter section
|
||||
# for higher distances use a linear vmin/vmax to cut out late/early regions with high noise
|
||||
if distance > 5.:
|
||||
pwleft = distance/vmax ################## TEST
|
||||
pwright = distance/vmin
|
||||
if pwright > cutwindow[1]:
|
||||
pwright = cutwindow[1]
|
||||
pickwin_used = (pwleft, pwright)
|
||||
|
||||
shot.setPickwindow(traceID, pickwin_used)
|
||||
shot.pickTraces(traceID, windowsize, folm, HosAic) # picker
|
||||
#shot.setManualPicks(traceID, picklist) # set manual picks if given (yet used on 2D only)
|
||||
|
||||
# ++ TEST: set and check SNR before adding to distance bin ############################
|
||||
shot.setSNR(traceID)
|
||||
#if shot.getSNR(traceID)[0] < snrthreshold:
|
||||
if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID):
|
||||
shot.removePick(traceID)
|
||||
# -- TEST: set and check SNR before adding to distance bin ############################
|
||||
|
||||
if shot.getPick(traceID) is not None:
|
||||
shot.setEarllatepick(traceID)
|
||||
|
||||
tpicksum += (datetime.now() - tstartpick); tpick = tpicksum/count
|
||||
tremain = (tpick * (len(survey.getShotDict()) - count))
|
||||
tend = datetime.now() + tremain
|
||||
print 'shot: %s, est. time to be finished is %s:%s:%s' % (shot.getShotname(), tend.hour, tend.minute, tend.second)
|
||||
|
||||
survey.saveSurvey(filename)
|
||||
print '\n--- Finished script ---'
|
||||
print 'Elapsed time:', datetime.now()-starttime
|
||||
@ -1,666 +0,0 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
import sys
|
||||
import numpy as np
|
||||
from scipy.interpolate import griddata
|
||||
|
||||
class SeisArray(object):
|
||||
'''
|
||||
Can be used to interpolate missing values of a receiver grid, if only support points were measured.
|
||||
Input file should contain in each line: ('traceID' 'receiverLineID' 'number of the geophone on recLine' 'X' 'Y' 'Z')
|
||||
|
||||
Can be used to generate a velocity grid file (vgrids.in) for FMTOMO with a topography adapting gradient.
|
||||
|
||||
Can be used to generate an interface file for FMTOMO (right now only interface.z used by grid3dg) for the topography.
|
||||
|
||||
Supports vtk output for sources and receivers.
|
||||
Note: Source and Receiver files for FMTOMO will be generated by the Survey object (because traveltimes will be added directly).
|
||||
'''
|
||||
def __init__(self, recfile):
|
||||
self._receiverlines = {}
|
||||
self._receiverCoords = {}
|
||||
self._measuredReceivers = {}
|
||||
self._measuredTopo = {}
|
||||
self._sourceLocs = {}
|
||||
self._geophoneNumbers = {}
|
||||
self._receiverlist = open(recfile, 'r').readlines()
|
||||
self._generateReceiverlines()
|
||||
self._setReceiverCoords()
|
||||
self._setGeophoneNumbers()
|
||||
|
||||
def _generateReceiverlines(self):
|
||||
'''
|
||||
Connects the traceIDs to the lineIDs
|
||||
for each receiverline in a dictionary.
|
||||
'''
|
||||
for receiver in self._receiverlist:
|
||||
traceID = int(receiver.split()[0])
|
||||
lineID = int(receiver.split()[1])
|
||||
if not lineID in self._receiverlines.keys():
|
||||
self._receiverlines[lineID] = []
|
||||
self._receiverlines[lineID].append(traceID)
|
||||
|
||||
def _setReceiverCoords(self):
|
||||
'''
|
||||
Fills the three x, y, z dictionaries with measured coordinates
|
||||
'''
|
||||
for line in self._getReceiverlist():
|
||||
traceID = int(line.split()[0])
|
||||
x = float(line.split()[3])
|
||||
y = float(line.split()[4])
|
||||
z = float(line.split()[5])
|
||||
self._receiverCoords[traceID] = (x, y, z)
|
||||
self._measuredReceivers[traceID] = (x, y, z)
|
||||
|
||||
def _setGeophoneNumbers(self):
|
||||
for line in self._getReceiverlist():
|
||||
traceID = int(line.split()[0])
|
||||
gphoneNum = float(line.split()[2])
|
||||
self._geophoneNumbers[traceID] = gphoneNum
|
||||
|
||||
def _getReceiverlines(self):
|
||||
return self._receiverlines
|
||||
|
||||
def _getReceiverlist(self):
|
||||
return self._receiverlist
|
||||
|
||||
def getReceiverCoordinates(self):
|
||||
return self._receiverCoords
|
||||
|
||||
def _getXreceiver(self, traceID):
|
||||
return self._receiverCoords[traceID][0]
|
||||
|
||||
def _getYreceiver(self, traceID):
|
||||
return self._receiverCoords[traceID][1]
|
||||
|
||||
def _getZreceiver(self, traceID):
|
||||
return self._receiverCoords[traceID][2]
|
||||
|
||||
def _getXshot(self, shotnumber):
|
||||
return self._sourceLocs[shotnumber][0]
|
||||
|
||||
def _getYshot(self, shotnumber):
|
||||
return self._sourceLocs[shotnumber][1]
|
||||
|
||||
def _getZshot(self, shotnumber):
|
||||
return self._sourceLocs[shotnumber][2]
|
||||
|
||||
def _getReceiverValue(self, traceID, coordinate):
|
||||
setCoordinate = {'X': self._getXreceiver,
|
||||
'Y': self._getYreceiver,
|
||||
'Z': self._getZreceiver}
|
||||
return setCoordinate[coordinate](traceID)
|
||||
|
||||
def _getGeophoneNumber(self, traceID):
|
||||
return self._geophoneNumbers[traceID]
|
||||
|
||||
def getMeasuredReceivers(self):
|
||||
return self._measuredReceivers
|
||||
|
||||
def getMeasuredTopo(self):
|
||||
return self._measuredTopo
|
||||
|
||||
def getSourceLocations(self):
|
||||
return self._sourceLocs
|
||||
|
||||
def _setXvalue(self, traceID, value):
|
||||
self._checkKey(traceID)
|
||||
self._receiverCoords[traceID][0] = value
|
||||
|
||||
def _setYvalue(self, traceID, value):
|
||||
self._checkKey(traceID)
|
||||
self._receiverCoords[traceID][1] = value
|
||||
|
||||
def _setZvalue(self, traceID, value):
|
||||
self._checkKey(traceID)
|
||||
self._receiverCoords[traceID][2] = value
|
||||
|
||||
def _setValue(self, traceID, coordinate, value):
|
||||
setCoordinate = {'X': self._setXvalue,
|
||||
'Y': self._setYvalue,
|
||||
'Z': self._setZvalue}
|
||||
setCoordinate[coordinate](traceID, value)
|
||||
|
||||
def _checkKey(self, traceID):
|
||||
if not traceID in self._receiverCoords.keys():
|
||||
self._receiverCoords[traceID] = [None, None, None]
|
||||
|
||||
def _checkTraceIDdirection(self, traceID1, traceID2):
|
||||
if traceID2 > traceID1:
|
||||
direction = +1
|
||||
return direction
|
||||
if traceID2 < traceID1:
|
||||
direction = -1
|
||||
return direction
|
||||
print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2)
|
||||
|
||||
def _checkCoordDirection(self, traceID1, traceID2, coordinate):
|
||||
'''
|
||||
Checks whether the interpolation direction is positive or negative
|
||||
'''
|
||||
if self._getReceiverValue(traceID1, coordinate) < self._getReceiverValue(traceID2, coordinate):
|
||||
direction = +1
|
||||
return direction
|
||||
if self._getReceiverValue(traceID1, coordinate) > self._getReceiverValue(traceID2, coordinate):
|
||||
direction = -1
|
||||
return direction
|
||||
print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2)
|
||||
|
||||
def _interpolateMeanDistances(self, traceID1, traceID2, coordinate):
|
||||
'''
|
||||
Returns the mean distance between two traceID's depending on the number of geophones in between
|
||||
'''
|
||||
num_spaces = abs(self._getGeophoneNumber(traceID1) - self._getGeophoneNumber(traceID2))
|
||||
mean_distance = abs(self._getReceiverValue(traceID1, coordinate) - self._getReceiverValue(traceID2, coordinate))/num_spaces
|
||||
return mean_distance
|
||||
|
||||
def interpolateValues(self, coordinate):
|
||||
'''
|
||||
Interpolates and sets all values (linear) for coordinate = 'X', 'Y' or 'Z'
|
||||
'''
|
||||
for lineID in self._getReceiverlines().keys():
|
||||
number_measured = len(self._getReceiverlines()[lineID])
|
||||
for index, traceID1 in enumerate(self._getReceiverlines()[lineID]):
|
||||
if index + 1 < number_measured:
|
||||
traceID2 = self._getReceiverlines()[lineID][index + 1]
|
||||
|
||||
traceID_dir = self._checkTraceIDdirection(traceID1, traceID2)
|
||||
traceID_interp = traceID1 + traceID_dir
|
||||
|
||||
coord_dir = self._checkCoordDirection(traceID1, traceID2, coordinate)
|
||||
mean_distance = self._interpolateMeanDistances(traceID1, traceID2, coordinate)
|
||||
|
||||
while (traceID_dir * traceID_interp) < (traceID_dir * traceID2):
|
||||
self._setValue(traceID_interp, coordinate,
|
||||
(self._getReceiverValue(traceID_interp - traceID_dir, coordinate)
|
||||
+ coord_dir * (mean_distance)))
|
||||
traceID_interp += traceID_dir
|
||||
|
||||
def addMeasuredTopographyPoints(self, filename):
|
||||
'''
|
||||
Use more measured points for a higher precision of height interpolation.
|
||||
Input file should contain in each line: ('point ID' 'X' 'Y' 'Z')
|
||||
'''
|
||||
topolist = open(filename, 'r').readlines()
|
||||
for line in topolist:
|
||||
line = line.split()
|
||||
pointID = int(line[0])
|
||||
x = float(line[1])
|
||||
y = float(line[2])
|
||||
z = float(line[3])
|
||||
self._measuredTopo[pointID] = (x, y, z)
|
||||
|
||||
def addSourceLocations(self, filename):
|
||||
'''
|
||||
Use source locations for a higher precision of height interpolation.
|
||||
Input file should contain in each line: ('point ID' 'X' 'Y' 'Z')
|
||||
|
||||
Source locations must be added before they can be written to vtk files.
|
||||
'''
|
||||
topolist = open(filename, 'r').readlines()
|
||||
for line in topolist:
|
||||
line = line.split()
|
||||
pointID = int(line[0])
|
||||
x = float(line[1])
|
||||
y = float(line[2])
|
||||
z = float(line[3])
|
||||
self._sourceLocs[pointID] = (x, y, z)
|
||||
|
||||
def interpZcoords4rec(self, method = 'linear'):
|
||||
'''
|
||||
Interpolates z values for all receivers.
|
||||
'''
|
||||
measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists()
|
||||
|
||||
for traceID in self.getReceiverCoordinates().keys():
|
||||
if type(self.getReceiverCoordinates()[traceID]) is not tuple:
|
||||
z = griddata((measured_x, measured_y), measured_z, (self._getXreceiver(traceID), self._getYreceiver(traceID)), method = method)
|
||||
self._setZvalue(traceID, float(z))
|
||||
|
||||
def _getAngle(self, distance):
|
||||
'''
|
||||
Function returns the angle on a Sphere of the radius R = 6371 [km] for a distance [km].
|
||||
'''
|
||||
PI = np.pi
|
||||
R = 6371.
|
||||
angle = distance * 180 / (PI * R)
|
||||
return angle
|
||||
|
||||
def _getDistance(self, angle):
|
||||
'''
|
||||
Function returns the distance [km] on a Sphere of the radius R = 6371 [km] for an angle.
|
||||
'''
|
||||
PI = np.pi
|
||||
R = 6371.
|
||||
distance = angle / 180 * (PI * R)
|
||||
return distance
|
||||
|
||||
def getMeasuredReceiverLists(self):
|
||||
'''
|
||||
Returns a list of all measured receivers known to SeisArray.
|
||||
'''
|
||||
x = []; y = []; z = []
|
||||
for traceID in self.getMeasuredReceivers().keys():
|
||||
x.append(self.getMeasuredReceivers()[traceID][0])
|
||||
y.append(self.getMeasuredReceivers()[traceID][1])
|
||||
z.append(self.getMeasuredReceivers()[traceID][2])
|
||||
return x, y, z
|
||||
|
||||
def getMeasuredTopoLists(self):
|
||||
'''
|
||||
Returns a list of all measured topography points known to the SeisArray.
|
||||
'''
|
||||
x = []; y = []; z = []
|
||||
for pointID in self.getMeasuredTopo().keys():
|
||||
x.append(self.getMeasuredTopo()[pointID][0])
|
||||
y.append(self.getMeasuredTopo()[pointID][1])
|
||||
z.append(self.getMeasuredTopo()[pointID][2])
|
||||
return x, y, z
|
||||
|
||||
def getSourceLocsLists(self):
|
||||
'''
|
||||
Returns a list of all measured source locations known to SeisArray.
|
||||
'''
|
||||
x = []; y = []; z = []
|
||||
for pointID in self.getSourceLocations().keys():
|
||||
x.append(self.getSourceLocations()[pointID][0])
|
||||
y.append(self.getSourceLocations()[pointID][1])
|
||||
z.append(self.getSourceLocations()[pointID][2])
|
||||
return x, y, z
|
||||
|
||||
def getAllMeasuredPointsLists(self):
|
||||
'''
|
||||
Returns a list of all measured points known to SeisArray.
|
||||
'''
|
||||
mtopo_x, mtopo_y, mtopo_z = self.getMeasuredTopoLists()
|
||||
msource_x, msource_y, msource_z = self.getSourceLocsLists()
|
||||
mrec_x, mrec_y, mrec_z = self.getMeasuredReceiverLists()
|
||||
|
||||
x = mtopo_x + mrec_x + msource_x
|
||||
y = mtopo_y + mrec_y + msource_y
|
||||
z = mtopo_z + mrec_z + msource_z
|
||||
return x, y, z
|
||||
|
||||
def getReceiverLists(self):
|
||||
'''
|
||||
Returns a list of all receivers (measured and interpolated).
|
||||
'''
|
||||
x = []; y =[]; z = []
|
||||
for traceID in self.getReceiverCoordinates().keys():
|
||||
x.append(self.getReceiverCoordinates()[traceID][0])
|
||||
y.append(self.getReceiverCoordinates()[traceID][1])
|
||||
z.append(self.getReceiverCoordinates()[traceID][2])
|
||||
return x, y, z
|
||||
|
||||
def _interpolateXY4rec(self):
|
||||
'''
|
||||
Interpolates the X and Y coordinates for all receivers.
|
||||
'''
|
||||
for coordinate in ('X', 'Y'):
|
||||
self.interpolateValues(coordinate)
|
||||
|
||||
def interpolateAll(self):
|
||||
self._interpolateXY4rec()
|
||||
self.interpZcoords4rec()
|
||||
|
||||
def interpolateTopography(self, nTheta, nPhi, thetaSN, phiWE, method = 'linear', filename = 'interface1.in'):
|
||||
'''
|
||||
Interpolate Z values on a regular grid with cushion nodes to use it as FMTOMO topography interface.
|
||||
Returns a surface in form of a list of points [[x1, y1, z1], [x2, y2, y2], ...] (cartesian).
|
||||
|
||||
:param: nTheta, number of points in theta (NS)
|
||||
type: integer
|
||||
|
||||
:param: nPhi, number of points in phi (WE)
|
||||
type: integer
|
||||
|
||||
:param: thetaSN (S, N) extensions of the model in degree
|
||||
type: tuple
|
||||
|
||||
:param: phiWE (W, E) extensions of the model in degree
|
||||
type: tuple
|
||||
'''
|
||||
|
||||
surface = []
|
||||
elevation = 0.25 # elevate topography so that no source lies above the surface
|
||||
|
||||
if filename is not None:
|
||||
outfile = open(filename, 'w')
|
||||
|
||||
print "Interpolating topography on regular grid with the dimensions:"
|
||||
print "nTheta = %s, nPhi = %s, thetaSN = %s, phiWE = %s"%(nTheta, nPhi, thetaSN, phiWE)
|
||||
print "method = %s, filename = %s" %(method, filename)
|
||||
|
||||
thetaS, thetaN = thetaSN
|
||||
phiW, phiE = phiWE
|
||||
|
||||
measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists()
|
||||
|
||||
# need to determine the delta to add two cushion nodes around the min/max values
|
||||
thetaDelta = (thetaN - thetaS) / (nTheta - 1)
|
||||
phiDelta = (phiE - phiW) / (nPhi - 1)
|
||||
|
||||
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes
|
||||
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes
|
||||
|
||||
nTotal = len(thetaGrid) * len(phiGrid); count = 0
|
||||
for theta in thetaGrid:
|
||||
for phi in phiGrid:
|
||||
xval = self._getDistance(phi)
|
||||
yval = self._getDistance(theta)
|
||||
z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = method)
|
||||
# in case the point lies outside, nan will be returned. Find nearest:
|
||||
if np.isnan(z) == True:
|
||||
z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = 'nearest')
|
||||
z = float(z)
|
||||
surface.append((xval, yval, z))
|
||||
count += 1
|
||||
progress = float(count) / float(nTotal) * 100
|
||||
self._update_progress(progress)
|
||||
|
||||
if filename is not None:
|
||||
outfile.writelines('%10s\n'%(z + elevation))
|
||||
|
||||
return surface
|
||||
|
||||
def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120,
|
||||
thetaSN = (-0.2, 1.2), phiWE = (-0.2, 1.2),
|
||||
Rbt = (-62.0, 6.0), vbot = 5.5, filename = 'vgrids.in',
|
||||
method = 'linear' ):
|
||||
'''
|
||||
Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s].
|
||||
|
||||
:param: nTheta, number of points in theta (NS)
|
||||
type: integer
|
||||
|
||||
:param: nPhi, number of points in phi (WE)
|
||||
type: integer
|
||||
|
||||
:param: nR, number of points in depth
|
||||
type: integer
|
||||
|
||||
:param: thetaSN (S, N) extensions of the model in degree
|
||||
type: tuple
|
||||
|
||||
:param: phiWE (W, E) extensions of the model in degree
|
||||
type: tuple
|
||||
|
||||
:param: Rbt (bot, top) extensions of the model in km
|
||||
type: tuple
|
||||
|
||||
:param: vbot, velocity at the bottom of the model
|
||||
type: real
|
||||
'''
|
||||
|
||||
def getRad(angle):
|
||||
PI = np.pi
|
||||
rad = angle / 180 * PI
|
||||
return rad
|
||||
|
||||
def getZmax(surface):
|
||||
z = []
|
||||
for point in surface:
|
||||
z.append(point[2])
|
||||
return max(z)
|
||||
|
||||
R = 6371
|
||||
vmin = 0.34
|
||||
decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
|
||||
outfile = open(filename, 'w')
|
||||
|
||||
thetaS, thetaN = thetaSN
|
||||
phiW, phiE = phiWE
|
||||
rbot = Rbt[0] + R
|
||||
rtop = Rbt[1] + R
|
||||
|
||||
# need to determine the delta to add two cushion nodes around the min/max values
|
||||
thetaDelta = abs(thetaN - thetaS) / (nTheta - 1)
|
||||
phiDelta = abs(phiE - phiW) / (nPhi - 1)
|
||||
rDelta = abs(rbot - rtop) / (nR - 1)
|
||||
|
||||
# create a regular grid including +2 cushion nodes in every direction
|
||||
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes
|
||||
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes
|
||||
rGrid = np.linspace(rbot - rDelta, rtop + rDelta, num = nR + 2) # +2 cushion nodes
|
||||
|
||||
nTotal = len(rGrid) * len(thetaGrid) * len(phiGrid)
|
||||
print "Total number of grid nodes: %s"%nTotal
|
||||
|
||||
# write header for velocity grid file (in RADIANS)
|
||||
outfile.writelines('%10s %10s \n' %(1, 1))
|
||||
outfile.writelines('%10s %10s %10s\n' %(nR + 2, nTheta + 2, nPhi + 2))
|
||||
outfile.writelines('%10s %10s %10s\n' %(rDelta, getRad(thetaDelta), getRad(phiDelta)))
|
||||
outfile.writelines('%10s %10s %10s\n' %(rbot - rDelta, getRad(thetaS - thetaDelta), getRad(phiW - phiDelta)))
|
||||
|
||||
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None)
|
||||
zmax = getZmax(surface)
|
||||
|
||||
print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(filename, method)
|
||||
print "nTheta = %s, nPhi = %s, nR = %s, thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt)
|
||||
count = 0
|
||||
for radius in rGrid:
|
||||
for theta in thetaGrid:
|
||||
for phi in phiGrid:
|
||||
xval = self._getDistance(phi)
|
||||
yval = self._getDistance(theta)
|
||||
for point in surface:
|
||||
if point[0] == xval and point[1] == yval:
|
||||
z = point[2]
|
||||
if radius > (R + z + 1):
|
||||
vel = 0.0
|
||||
# elif radius > (R + z - 15): ########### TESTING
|
||||
# vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * 1.0 + vmin ##########################
|
||||
else:
|
||||
vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * vbot + vmin ##########################
|
||||
count += 1
|
||||
outfile.writelines('%10s %10s\n'%(vel, decm))
|
||||
|
||||
progress = float(count) / float(nTotal) * 100
|
||||
self._update_progress(progress)
|
||||
|
||||
outfile.close()
|
||||
|
||||
def exportAll(self, filename = 'interpolated_receivers.out'):
|
||||
recfile_out = open(filename, 'w')
|
||||
count = 0
|
||||
for traceID in self.getReceiverCoordinates().keys():
|
||||
count += 1
|
||||
x, y, z = self.getReceiverCoordinates()[traceID]
|
||||
recfile_out.writelines('%5s %15s %15s %15s\n' %(traceID, x, y, z))
|
||||
print "Exported coordinates for %s traces to file > %s" %(count, filename)
|
||||
recfile_out.close()
|
||||
|
||||
def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True):
|
||||
import matplotlib.pyplot as plt
|
||||
plt.interactive(True)
|
||||
plt.figure()
|
||||
xmt, ymt, zmt = self.getMeasuredTopoLists()
|
||||
xsc, ysc, zsc = self.getSourceLocsLists()
|
||||
xmr, ymr, zmr = self.getMeasuredReceiverLists()
|
||||
xrc, yrc, zrc = self.getReceiverLists()
|
||||
|
||||
plt.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers')
|
||||
plt.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations')
|
||||
|
||||
if plot_topo == True:
|
||||
plt.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points')
|
||||
if highlight_measured == True:
|
||||
plt.plot(xmr, ymr, 'ro', label = 'measured receivers')
|
||||
|
||||
plt.xlabel('X [m]')
|
||||
plt.ylabel('Y [m]')
|
||||
plt.legend()
|
||||
if annotations == True:
|
||||
for traceID in self.getReceiverCoordinates().keys():
|
||||
plt.annotate(str(traceID), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small')
|
||||
|
||||
def plotArray3D(self, ax = None):
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
plt.interactive(True)
|
||||
|
||||
if ax == None:
|
||||
fig = plt.figure()
|
||||
ax = plt.axes(projection = '3d')
|
||||
|
||||
xmt, ymt, zmt = self.getMeasuredTopoLists()
|
||||
xmr, ymr, zmr = self.getMeasuredReceiverLists()
|
||||
xin, yin, zin = self.getReceiverLists()
|
||||
|
||||
ax.plot(xmt, ymt, zmt, 'b*', markersize = 10, label = 'measured topo points')
|
||||
ax.plot(xin, yin, zin, 'k.', markersize = 10, label = 'interpolated receivers')
|
||||
ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers')
|
||||
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation')
|
||||
ax.legend()
|
||||
|
||||
return ax
|
||||
|
||||
|
||||
def plotSurface3D(self, ax = None, step = 0.5, method = 'linear'):
|
||||
from matplotlib import cm
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
plt.interactive(True)
|
||||
|
||||
if ax == None:
|
||||
fig = plt.figure()
|
||||
ax = plt.axes(projection = '3d')
|
||||
|
||||
xmt, ymt, zmt = self.getMeasuredTopoLists()
|
||||
xmr, ymr, zmr = self.getMeasuredReceiverLists()
|
||||
|
||||
x = xmt + xmr
|
||||
y = ymt + ymr
|
||||
z = zmt + zmr
|
||||
|
||||
xaxis = np.arange(min(x)+1, max(x), step)
|
||||
yaxis = np.arange(min(y)+1, max(y), step)
|
||||
|
||||
xgrid, ygrid = np.meshgrid(xaxis, yaxis)
|
||||
|
||||
zgrid = griddata((x, y), z, (xgrid, ygrid), method = method)
|
||||
|
||||
ax.plot_surface(xgrid, ygrid, zgrid, linewidth = 0, cmap = cm.jet, vmin = min(z), vmax = max(z))
|
||||
|
||||
ax.set_zlim(-(max(x) - min(x)/2),(max(x) - min(x)/2))
|
||||
ax.set_aspect('equal')
|
||||
|
||||
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation')
|
||||
ax.legend()
|
||||
|
||||
return ax
|
||||
|
||||
def _update_progress(self, progress):
|
||||
sys.stdout.write("%d%% done \r" % (progress) )
|
||||
sys.stdout.flush()
|
||||
|
||||
def receivers2VTK(self, filename = 'receivers.vtk'):
|
||||
'''
|
||||
Generates vtk files from all receivers of the SeisArray object.
|
||||
'''
|
||||
outfile = open(filename, 'w')
|
||||
traceIDs = []
|
||||
|
||||
for traceID in self.getReceiverCoordinates():
|
||||
traceIDs.append(traceID)
|
||||
|
||||
nPoints = len(traceIDs)
|
||||
|
||||
# write header
|
||||
print("Writing header for VTK file...")
|
||||
outfile.writelines('# vtk DataFile Version 3.1\n')
|
||||
outfile.writelines('Receivers with traceIDs\n')
|
||||
outfile.writelines('ASCII\n')
|
||||
outfile.writelines('DATASET POLYDATA\n')
|
||||
outfile.writelines('POINTS %15d float\n' %(nPoints))
|
||||
|
||||
# write coordinates
|
||||
print("Writing coordinates to VTK file...")
|
||||
for traceID in traceIDs:
|
||||
x = self._getXreceiver(traceID)
|
||||
y = self._getYreceiver(traceID)
|
||||
z = self._getZreceiver(traceID)
|
||||
|
||||
outfile.writelines('%10f %10f %10f \n' %(x, y, z))
|
||||
|
||||
outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
|
||||
|
||||
# write indices
|
||||
print("Writing indices to VTK file...")
|
||||
for index in range(nPoints):
|
||||
outfile.writelines('%10d %10d\n' %(1, index))
|
||||
|
||||
outfile.writelines('POINT_DATA %15d\n' %(nPoints))
|
||||
outfile.writelines('SCALARS traceIDs int %d\n' %(1))
|
||||
outfile.writelines('LOOKUP_TABLE default\n')
|
||||
|
||||
# write traceIDs
|
||||
print("Writing traceIDs to VTK file...")
|
||||
for traceID in traceIDs:
|
||||
outfile.writelines('%10d\n' %traceID)
|
||||
|
||||
outfile.close()
|
||||
print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename))
|
||||
return
|
||||
|
||||
def sources2VTK(self, filename = 'sources.vtk'):
|
||||
'''
|
||||
Generates vtk-files for all source locations in the SeisArray object.
|
||||
'''
|
||||
outfile = open(filename, 'w')
|
||||
shotnumbers = []
|
||||
|
||||
for shotnumber in self.getSourceLocations():
|
||||
shotnumbers.append(shotnumber)
|
||||
|
||||
nPoints = len(shotnumbers)
|
||||
|
||||
# write header
|
||||
print("Writing header for VTK file...")
|
||||
outfile.writelines('# vtk DataFile Version 3.1\n')
|
||||
outfile.writelines('Shots with shotnumbers\n')
|
||||
outfile.writelines('ASCII\n')
|
||||
outfile.writelines('DATASET POLYDATA\n')
|
||||
outfile.writelines('POINTS %15d float\n' %(nPoints))
|
||||
|
||||
# write coordinates
|
||||
print("Writing coordinates to VTK file...")
|
||||
for shotnumber in shotnumbers:
|
||||
x = self._getXshot(shotnumber)
|
||||
y = self._getYshot(shotnumber)
|
||||
z = self._getZshot(shotnumber)
|
||||
|
||||
outfile.writelines('%10f %10f %10f \n' %(x, y, z))
|
||||
|
||||
outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
|
||||
|
||||
# write indices
|
||||
print("Writing indices to VTK file...")
|
||||
for index in range(nPoints):
|
||||
outfile.writelines('%10d %10d\n' %(1, index))
|
||||
|
||||
outfile.writelines('POINT_DATA %15d\n' %(nPoints))
|
||||
outfile.writelines('SCALARS shotnumbers int %d\n' %(1))
|
||||
outfile.writelines('LOOKUP_TABLE default\n')
|
||||
|
||||
# write shotnumber
|
||||
print("Writing shotnumbers to VTK file...")
|
||||
for shotnumber in shotnumbers:
|
||||
outfile.writelines('%10d\n' %shotnumber)
|
||||
|
||||
outfile.close()
|
||||
print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename))
|
||||
return
|
||||
|
||||
|
||||
def saveSeisArray(self, filename = 'seisArray.pickle'):
|
||||
import cPickle
|
||||
outfile = open(filename, 'wb')
|
||||
|
||||
cPickle.dump(self, outfile, -1)
|
||||
print('saved SeisArray to file %s'%(filename))
|
||||
|
||||
@staticmethod
|
||||
def from_pickle(filename):
|
||||
import cPickle
|
||||
infile = open(filename, 'rb')
|
||||
return cPickle.load(infile)
|
||||
Loading…
Reference in New Issue
Block a user