1051 lines
39 KiB
Python
1051 lines
39 KiB
Python
# -*- 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.recfile = recfile
|
|
self._receiverlines = {}
|
|
self._receiverCoords = {}
|
|
self._measuredReceivers = {}
|
|
self._measuredTopo = {}
|
|
self._sourceLocs = {}
|
|
self._geophoneNumbers = {}
|
|
self._receiverlist = open(self.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
|
|
err_msg = "Same Value for traceID1 = %s and traceID2 = %s" % (traceID1, traceID2)
|
|
raise RuntimeError(err_msg)
|
|
|
|
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
|
|
err_msg = "Same Value for traceID1 = %s and traceID2 = %s" % (traceID1, traceID2)
|
|
raise RuntimeError(err_msg)
|
|
|
|
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()
|
|
print('Interpolated receiver locations.')
|
|
|
|
def interpolateTopography(self, nTheta, nPhi, thetaSN, phiWE, elevation=0.25, method='linear'):
|
|
'''
|
|
Interpolate Z values on a regular grid with cushion nodes e.g. 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
|
|
|
|
:param: elevation, default: 0.25 (elevate topography so that no source lies above the surface)
|
|
type: float
|
|
'''
|
|
return self.interpolateOnRegularGrid(nTheta, nPhi, thetaSN, phiWE, elevation, method)
|
|
|
|
def interpolateOnRegularGrid(self, nTheta, nPhi, thetaSN, phiWE, elevation, method='linear'):
|
|
'''
|
|
Interpolate Z values on a regular grid with cushion nodes e.g. 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
|
|
|
|
:param: elevation, default: 0.25 (elevate topography so that no source lies above the surface)
|
|
type: float
|
|
'''
|
|
|
|
surface = []
|
|
|
|
print "Interpolating interface on regular grid with the dimensions:"
|
|
print "nTheta = %s, nPhi = %s, thetaSN = %s, phiWE = %s" % (nTheta, nPhi, thetaSN, phiWE)
|
|
print "method = %s, elevation = %s" % (method, elevation)
|
|
|
|
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
|
|
deltaTheta = (thetaN - thetaS) / (nTheta - 1)
|
|
deltaPhi = (phiE - phiW) / (nPhi - 1)
|
|
|
|
thetaGrid = np.linspace(thetaS - deltaTheta, thetaN + deltaTheta, num=nTheta + 2) # +2 cushion nodes
|
|
phiGrid = np.linspace(phiW - deltaPhi, phiE + deltaPhi, 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) + elevation
|
|
surface.append((xval, yval, z))
|
|
count += 1
|
|
progress = float(count) / float(nTotal) * 100
|
|
self._update_progress(progress)
|
|
|
|
return surface
|
|
|
|
def generateFMTOMOinputFromArray(self, nPointsPropgrid, nPointsInvgrid,
|
|
zBotTop, cushionfactor, interpolationMethod='linear',
|
|
customgrid='mygrid.in', writeVTK=True):
|
|
'''
|
|
Generate FMTOMO input files from the SeisArray dimensions.
|
|
Generates: vgridsref.in, interfacesref.in, propgrid.in
|
|
|
|
:param: nPointsPropgrid, number of points in each direction of the propagation grid (z, y, x)
|
|
:type: tuple
|
|
|
|
:param: nPointsInvgrid, number of points in each direction of the inversion grid (z, y, x)
|
|
:type: tuple
|
|
|
|
:param: zBotTop, (bottom, top) dimensions of the model
|
|
:type: tuple
|
|
|
|
:param: cushionfactor, adds cushioning around the model (0.1 = 10%)
|
|
:type: float
|
|
'''
|
|
|
|
nRP, nThetaP, nPhiP = nPointsPropgrid
|
|
nRI, nThetaI, nPhiI = nPointsInvgrid
|
|
|
|
print('\n------------------------------------------------------------')
|
|
print('Automatically generating input for FMTOMO from array size.')
|
|
print('Propgrid: nR = %s, nTheta = %s, nPhi = %s' % (nRP, nThetaP, nPhiP))
|
|
print('Interpolation Grid and Interfaces Grid: nR = %s, nTheta = %s, nPhi = %s' % (nRI, nThetaI, nPhiI))
|
|
print('Bottom and Top of model: (%s, %s)' % (zBotTop[0], zBotTop[1]))
|
|
print('Method: %s, customgrid = %s' % (interpolationMethod, customgrid))
|
|
print('------------------------------------------------------------')
|
|
|
|
def getZmin(surface):
|
|
z = []
|
|
for point in surface:
|
|
z.append(point[2])
|
|
return min(z)
|
|
|
|
self.generatePropgrid(nThetaP, nPhiP, nRP, zBotTop, cushionfactor=cushionfactor,
|
|
cushionpropgrid=0.05)
|
|
surface = self.generateVgrid(nThetaI, nPhiI, nRI, zBotTop, method=interpolationMethod,
|
|
cushionfactor=cushionfactor, infilename=customgrid,
|
|
returnTopo=True)
|
|
|
|
depthmax = abs(zBotTop[0] - getZmin(surface)) - 1.0 # cushioning for the bottom interface
|
|
|
|
interf1, interf2 = self.generateInterfaces(nThetaI, nPhiI, depthmax, cushionfactor=cushionfactor,
|
|
returnInterfaces=True, method=interpolationMethod)
|
|
|
|
if writeVTK == True:
|
|
from pylot.core.active import fmtomoUtils
|
|
self.surface2VTK(interf1, filename='interface1.vtk')
|
|
self.surface2VTK(interf2, filename='interface2.vtk')
|
|
self.receivers2VTK()
|
|
self.sources2VTK()
|
|
fmtomoUtils.vgrids2VTK()
|
|
|
|
def generateReceiversIn(self, outfilename='receivers.in'):
|
|
outfile = open(outfilename, 'w')
|
|
|
|
recx, recy, recz = self.getReceiverLists()
|
|
nsrc = len(self.getSourceLocations())
|
|
outfile.write('%s\n' % (len(zip(recx, recy, recz)) * nsrc))
|
|
|
|
for index in range(nsrc):
|
|
for point in zip(recx, recy, recz):
|
|
rx, ry, rz = point
|
|
rad = - rz
|
|
lat = self._getAngle(ry)
|
|
lon = self._getAngle(rx)
|
|
outfile.write('%15s %15s %15s\n' % (rad, lat, lon))
|
|
outfile.write('%15s\n' % (1))
|
|
outfile.write('%15s\n' % (index + 1))
|
|
outfile.write('%15s\n' % (1))
|
|
|
|
outfile.close()
|
|
|
|
def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor=0.1,
|
|
outfilename='interfacesref.in', method='linear',
|
|
returnInterfaces=False):
|
|
'''
|
|
Create an interfacesref.in file for FMTOMO from the SeisArray boundaries.
|
|
:param: nTheta, number of points in Theta
|
|
type: int
|
|
|
|
:param: nPhi, number of points in Phi
|
|
type: int
|
|
|
|
:param: depthmax, maximum depth of the model (below topography)
|
|
type: float
|
|
|
|
:param: cushionfactor, add some extra space to the model (default: 0.1 = 10%)
|
|
type: float
|
|
'''
|
|
|
|
print('\n------------------------------------------------------------')
|
|
print('Generating interfaces...')
|
|
nInterfaces = 2
|
|
|
|
# generate dimensions of the grid from array
|
|
thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor)
|
|
|
|
thetaS, thetaN = thetaSN
|
|
phiW, phiE = phiWE
|
|
R = 6371.
|
|
|
|
outfile = open(outfilename, 'w')
|
|
|
|
# determine the deltas
|
|
deltaTheta = abs(thetaN - thetaS) / float((nTheta - 1))
|
|
deltaPhi = abs(phiE - phiW) / float((nPhi - 1))
|
|
|
|
# write header for interfaces grid file (in RADIANS)
|
|
outfile.write('%10s\n' % (nInterfaces))
|
|
outfile.write('%10s %10s\n' % (nTheta + 2, nPhi + 2)) # +2 cushion nodes
|
|
outfile.write('%10s %10s\n' % (np.deg2rad(deltaTheta), np.deg2rad(deltaPhi)))
|
|
outfile.write('%10s %10s\n' % (np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi)))
|
|
|
|
interface1 = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method=method)
|
|
interface2 = self.interpolateOnRegularGrid(nTheta, nPhi, thetaSN, phiWE, -depthmax, method=method)
|
|
|
|
for point in interface1:
|
|
z = point[2]
|
|
outfile.write('%10s\n' % (z + R))
|
|
|
|
outfile.write('\n')
|
|
for point in interface2:
|
|
z = point[2]
|
|
outfile.write('%10s\n' % (z + R))
|
|
|
|
outfile.close()
|
|
|
|
if returnInterfaces == True:
|
|
return interface1, interface2
|
|
|
|
print('Finished generating interfaces.')
|
|
print('------------------------------------------------------------')
|
|
|
|
def getThetaPhiFromArray(self, cushionfactor=0.1):
|
|
'''
|
|
Determine and returns PhiWE (tuple: (West, East)) and thetaSN (tuple (South, North)) from the SeisArray boundaries.
|
|
|
|
:param: cushionfactor, add some extra space to the model (default: 0.1 = 10%)
|
|
type: float
|
|
'''
|
|
x, y, _ = self.getAllMeasuredPointsLists()
|
|
phi_min, phi_max = (self._getAngle(min(x)), self._getAngle(max(x)))
|
|
theta_min, theta_max = (self._getAngle(min(y)), self._getAngle(max(y)))
|
|
cushionPhi = abs(phi_max - phi_min) * cushionfactor
|
|
cushionTheta = abs(theta_max - theta_min) * cushionfactor
|
|
phiWE = (phi_min - cushionPhi, phi_max + cushionPhi)
|
|
thetaSN = (theta_min - cushionTheta, theta_max + cushionTheta)
|
|
return thetaSN, phiWE
|
|
|
|
def generatePropgrid(self, nTheta, nPhi, nR, Rbt, cushionfactor, cushionpropgrid=0.05,
|
|
refinement=(5, 5), outfilename='propgrid.in'):
|
|
'''
|
|
Create a propergation grid file for FMTOMO using SeisArray boundaries
|
|
|
|
:param: nTheta, number of points in Theta
|
|
type: int
|
|
|
|
:param: nPhi, number of points in Phi
|
|
type: int
|
|
|
|
:param: nR, number of points in R
|
|
type: int
|
|
|
|
:param: Rbt (bot, top) extensions of the model in m
|
|
type: tuple
|
|
|
|
:param: cushionfactor, add some extra space to the model (default: 0.1 = 10%)
|
|
type: float
|
|
|
|
:param: cushionpropogrid, cushionfactor for the propagationgrid (cushion direction
|
|
opposing to vgrids cushionfactor)
|
|
type: float
|
|
|
|
:param: refinement, (refinement factor, number of local cells for refinement) used by FMTOMO
|
|
type: tuple
|
|
'''
|
|
outfile = open(outfilename, 'w')
|
|
|
|
print('\n------------------------------------------------------------')
|
|
print('Generating Propagation Grid for nTheta = %s, nPhi'
|
|
' = %s, nR = %s and a cushioning of %s'
|
|
% (nTheta, nPhi, nR, cushionpropgrid))
|
|
print('Bottom of the grid: %s, top of the grid %s'
|
|
% (Rbt[0], Rbt[1]))
|
|
|
|
thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor)
|
|
|
|
thetaS = thetaSN[0] + cushionpropgrid
|
|
thetaN = thetaSN[1] - cushionpropgrid
|
|
phiW = phiWE[0] + cushionpropgrid
|
|
phiE = phiWE[1] - cushionpropgrid
|
|
rbot = Rbt[0]
|
|
rtop = Rbt[1]
|
|
|
|
deltaTheta = abs(thetaN - thetaS) / float(nTheta - 1)
|
|
deltaPhi = abs(phiE - phiW) / float(nPhi - 1)
|
|
deltaR = abs(rbot - rtop) / float(nR - 1)
|
|
|
|
outfile.write('%10s %10s %10s\n' % (nR, nTheta, nPhi))
|
|
outfile.write('%10s %10s %10s\n' % (deltaR, deltaTheta, deltaPhi))
|
|
outfile.write('%10s %10s %10s\n' % (rtop, thetaS, phiW))
|
|
outfile.write('%10s %10s\n' % refinement)
|
|
|
|
outfile.close()
|
|
|
|
print('Created Propagation Grid and saved it to %s' % outfilename)
|
|
print('------------------------------------------------------------')
|
|
|
|
def generateVgrid(self, nTheta, nPhi, nR, Rbt, thetaSN=None,
|
|
phiWE=None, cushionfactor=0.1,
|
|
outfilename='vgridsref.in', method='linear',
|
|
infilename='mygrid.in', returnTopo=False):
|
|
'''
|
|
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 m
|
|
type: tuple
|
|
|
|
:param: vbot, velocity at the bottom of the model
|
|
type: real
|
|
|
|
:param: method, interpolation method for topography
|
|
type: str
|
|
'''
|
|
print('\n------------------------------------------------------------')
|
|
print('generateVgrid: Starting...')
|
|
|
|
# def getRad(angle):
|
|
# PI = np.pi
|
|
# rad = angle / 180 * PI
|
|
# return rad
|
|
|
|
def readMygridNlayers(filename):
|
|
infile = open(filename, 'r')
|
|
nlayers = len(infile.readlines()) / 2
|
|
infile.close()
|
|
|
|
return nlayers
|
|
|
|
def readMygrid(filename):
|
|
ztop = [];
|
|
zbot = [];
|
|
vtop = [];
|
|
vbot = []
|
|
infile = open(filename, 'r')
|
|
nlayers = readMygridNlayers(filename)
|
|
|
|
print('\nreadMygrid: Reading file %s.' % filename)
|
|
for index in range(nlayers):
|
|
line1 = infile.readline()
|
|
line2 = infile.readline()
|
|
ztop.append(float(line1.split()[0]))
|
|
vtop.append(float(line1.split()[1]))
|
|
zbot.append(float(line2.split()[0]))
|
|
vbot.append(float(line2.split()[1]))
|
|
print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]'
|
|
'\n[Bot: v = %s [km/s], z = %s [m]]'
|
|
% (index + 1, vtop[index], ztop[index],
|
|
vbot[index], zbot[index]))
|
|
|
|
if not ztop[0] == 0:
|
|
print('ERROR: there must be a velocity set for z = 0 in the file %s' % filename)
|
|
print('e.g.:\n0 0.33\n-5 1.0\netc.')
|
|
|
|
infile.close()
|
|
return ztop, zbot, vtop, vbot
|
|
|
|
R = 6371.
|
|
vmin = 0.34
|
|
decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
|
|
outfile = open(outfilename, 'w')
|
|
|
|
# generate dimensions of the grid from array
|
|
if thetaSN is None and phiWE is None:
|
|
thetaSN, phiWE = self.getThetaPhiFromArray(cushionfactor)
|
|
|
|
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
|
|
deltaTheta = abs(thetaN - thetaS) / float((nTheta - 1))
|
|
deltaPhi = abs(phiE - phiW) / float((nPhi - 1))
|
|
deltaR = abs(rbot - rtop) / float((nR - 1))
|
|
|
|
# create a regular grid including +2 cushion nodes in every direction
|
|
thetaGrid = np.linspace(thetaS - deltaTheta, thetaN + deltaTheta, num=nTheta + 2) # +2 cushion nodes
|
|
phiGrid = np.linspace(phiW - deltaPhi, phiE + deltaPhi, num=nPhi + 2) # +2 cushion nodes
|
|
rGrid = np.linspace(rbot - deltaR, rtop + deltaR, 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.write('%10s %10s \n' % (1, 1))
|
|
outfile.write('%10s %10s %10s\n' % (nR + 2, nTheta + 2, nPhi + 2))
|
|
outfile.write('%10s %10s %10s\n' % (deltaR, np.deg2rad(deltaTheta), np.deg2rad(deltaPhi)))
|
|
outfile.write(
|
|
'%10s %10s %10s\n' % (rbot - deltaR, np.deg2rad(thetaS - deltaTheta), np.deg2rad(phiW - deltaPhi)))
|
|
|
|
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method=method)
|
|
|
|
nlayers = readMygridNlayers(infilename)
|
|
ztop, zbot, vtop, vbot = readMygrid(infilename)
|
|
|
|
print("\nGenerating velocity grid for FMTOMO. "
|
|
"Output filename = %s, interpolation method = %s" % (outfilename, 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:
|
|
topo = point[2]
|
|
break
|
|
depth = -(R + topo - radius)
|
|
if depth > 1:
|
|
vel = 0.0
|
|
elif 1 >= depth > 0: # cushioning around topography
|
|
vel = vtop[0]
|
|
else:
|
|
for index in range(nlayers):
|
|
if (ztop[index]) >= depth > (zbot[index]):
|
|
vel = (depth - ztop[index]) / (zbot[index] - ztop[index]) * (
|
|
vbot[index] - vtop[index]) + vtop[index]
|
|
break
|
|
if not (ztop[index]) >= depth > (zbot[index]):
|
|
err_msg = 'ERROR in grid inputfile, could not find velocity'
|
|
'for a z-value of %s in the inputfile' % (depth - topo)
|
|
raise ValueError(err_msg)
|
|
count += 1
|
|
if not vel >= 0:
|
|
err_msg = 'vel < 0; z, topo, zbot, vbot, vtop:',
|
|
depth, topo, zbot[index], vbot[index], vtop[index]
|
|
raise ValueError(err_msg)
|
|
|
|
outfile.write('%10s %10s\n' % (vel, decm))
|
|
|
|
progress = float(count) / float(nTotal) * 100
|
|
self._update_progress(progress)
|
|
|
|
print('\nWrote %d points to file %s for %d layers' % (count, outfilename, nlayers))
|
|
print('------------------------------------------------------------')
|
|
outfile.close()
|
|
|
|
if returnTopo == True:
|
|
return surface
|
|
|
|
def exportAll(self, filename='interpolated_receivers.out'):
|
|
'''
|
|
Exports all receivers to an input file for ActiveSeismoPick3D.
|
|
'''
|
|
recfile_out = open(filename, 'w')
|
|
count = 0
|
|
for traceID in self.getReceiverCoordinates().keys():
|
|
count += 1
|
|
x, y, z = self.getReceiverCoordinates()[traceID]
|
|
recfile_out.write('%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, ax = None, plot_topo=False, highlight_measured=False, annotations=True, pointsize=10):
|
|
import matplotlib.pyplot as plt
|
|
if ax == None:
|
|
plt.interactive(True)
|
|
fig = plt.figure()
|
|
ax = plt.axes()
|
|
xmt, ymt, zmt = self.getMeasuredTopoLists()
|
|
xsc, ysc, zsc = self.getSourceLocsLists()
|
|
xmr, ymr, zmr = self.getMeasuredReceiverLists()
|
|
xrc, yrc, zrc = self.getReceiverLists()
|
|
|
|
if len(xrc) > 0:
|
|
ax.plot(xrc, yrc, 'k.', markersize=pointsize, label='all receivers')
|
|
if len(xsc) > 0:
|
|
ax.plot(xsc, ysc, 'b*', markersize=pointsize, label='shot locations')
|
|
|
|
if plot_topo == True:
|
|
ax.plot(xmt, ymt, 'b.', markersize=pointsize, label='measured topo points')
|
|
if highlight_measured == True:
|
|
ax.plot(xmr, ymr, 'r.', markersize=pointsize, label='measured receivers')
|
|
|
|
ax.text(0.5, 1.05,'2D plot of seismic array\n %s'%self.recfile,
|
|
horizontalalignment='center', verticalalignment='center',
|
|
transform=ax.transAxes)
|
|
#plt.title('2D plot of seismic array %s' % self.recfile)
|
|
ax.set_xlabel('X [m]')
|
|
ax.set_ylabel('Y [m]')
|
|
ax.set_aspect('equal')
|
|
ax.legend()
|
|
if annotations == True:
|
|
for traceID in self.getReceiverCoordinates().keys():
|
|
ax.annotate((' ' + str(traceID)), xy=(self._getXreceiver(traceID), self._getYreceiver(traceID)),
|
|
fontsize='x-small', color='k')
|
|
for shotnumber in self.getSourceLocations().keys():
|
|
ax.annotate((' ' + str(shotnumber)), xy=(self._getXshot(shotnumber), self._getYshot(shotnumber)),
|
|
fontsize='x-small', color='b')
|
|
|
|
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()
|
|
xrc, yrc, zrc = self.getReceiverLists()
|
|
xsc, ysc, zsc = self.getSourceLocsLists()
|
|
|
|
plt.title('3D plot of seismic array.')
|
|
if len(xmt) > 0:
|
|
ax.plot(xmt, ymt, zmt, 'b.', markersize=10, label='measured topo points')
|
|
if len(xrc) > 0:
|
|
ax.plot(xrc, yrc, zrc, 'k.', markersize=10, label='all receivers')
|
|
if len(xmr) > 0:
|
|
ax.plot(xmr, ymr, zmr, 'ro', label='measured receivers')
|
|
if len(xsc) > 0:
|
|
ax.plot(xsc, ysc, zsc, 'b*', label='shot locations')
|
|
ax.set_xlabel('X [m]');
|
|
ax.set_ylabel('Y [m]');
|
|
ax.set_zlabel('Z [m]')
|
|
ax.legend()
|
|
|
|
return ax
|
|
|
|
def plotSurface3D(self, ax=None, step=0.5, method='linear', exag=False):
|
|
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)
|
|
|
|
surf = ax.plot_surface(xgrid, ygrid, zgrid, linewidth=0, cmap=cm.jet, vmin=min(z), vmax=max(z))
|
|
cbar = plt.colorbar(surf)
|
|
cbar.set_label('Elevation [m]')
|
|
|
|
if exag == False:
|
|
ax.set_zlim(-(max(x) - min(x) / 2), (max(x) - min(x) / 2))
|
|
ax.set_aspect('equal')
|
|
|
|
ax.set_xlabel('X [m]');
|
|
ax.set_ylabel('Y [m]');
|
|
ax.set_zlabel('Z [m]')
|
|
ax.legend()
|
|
|
|
return ax
|
|
|
|
def _update_progress(self, progress):
|
|
sys.stdout.write("%d%% done \r" % (progress))
|
|
sys.stdout.flush()
|
|
|
|
def surface2VTK(self, surface, filename='surface.vtk'):
|
|
'''
|
|
Generates a vtk file from all points of a surface as generated by interpolateTopography.
|
|
'''
|
|
outfile = open(filename, 'w')
|
|
|
|
nPoints = len(surface)
|
|
|
|
# write header
|
|
print("Writing header for VTK file...")
|
|
outfile.write('# vtk DataFile Version 3.1\n')
|
|
outfile.write('Surface Points\n')
|
|
outfile.write('ASCII\n')
|
|
outfile.write('DATASET POLYDATA\n')
|
|
outfile.write('POINTS %15d float\n' % (nPoints))
|
|
|
|
# write coordinates
|
|
print("Writing coordinates to VTK file...")
|
|
for point in surface:
|
|
x = point[0]
|
|
y = point[1]
|
|
z = point[2]
|
|
|
|
outfile.write('%10f %10f %10f \n' % (x, y, z))
|
|
|
|
outfile.write('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints))
|
|
|
|
# write indices
|
|
print("Writing indices to VTK file...")
|
|
for index in range(nPoints):
|
|
outfile.write('%10d %10d\n' % (1, index))
|
|
|
|
# outfile.write('POINT_DATA %15d\n' %(nPoints))
|
|
# outfile.write('SCALARS traceIDs int %d\n' %(1))
|
|
# outfile.write('LOOKUP_TABLE default\n')
|
|
|
|
# # write traceIDs
|
|
# print("Writing traceIDs to VTK file...")
|
|
# for traceID in traceIDs:
|
|
# outfile.write('%10d\n' %traceID)
|
|
|
|
outfile.close()
|
|
print("Wrote %d points to file: %s" % (nPoints, filename))
|
|
return
|
|
|
|
def receivers2VTK(self, filename='receivers.vtk'):
|
|
'''
|
|
Generates a vtk file 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.write('# vtk DataFile Version 3.1\n')
|
|
outfile.write('Receivers with traceIDs\n')
|
|
outfile.write('ASCII\n')
|
|
outfile.write('DATASET POLYDATA\n')
|
|
outfile.write('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.write('%10f %10f %10f \n' % (x, y, z))
|
|
|
|
outfile.write('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints))
|
|
|
|
# write indices
|
|
print("Writing indices to VTK file...")
|
|
for index in range(nPoints):
|
|
outfile.write('%10d %10d\n' % (1, index))
|
|
|
|
outfile.write('POINT_DATA %15d\n' % (nPoints))
|
|
outfile.write('SCALARS traceIDs int %d\n' % (1))
|
|
outfile.write('LOOKUP_TABLE default\n')
|
|
|
|
# write traceIDs
|
|
print("Writing traceIDs to VTK file...")
|
|
for traceID in traceIDs:
|
|
outfile.write('%10d\n' % traceID)
|
|
|
|
outfile.close()
|
|
print("Wrote %d receiver for to file: %s" % (nPoints, filename))
|
|
return
|
|
|
|
def sources2VTK(self, filename='sources.vtk'):
|
|
'''
|
|
Generates a vtk-file 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.write('# vtk DataFile Version 3.1\n')
|
|
outfile.write('Shots with shotnumbers\n')
|
|
outfile.write('ASCII\n')
|
|
outfile.write('DATASET POLYDATA\n')
|
|
outfile.write('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.write('%10f %10f %10f \n' % (x, y, z))
|
|
|
|
outfile.write('VERTICES %15d %15d\n' % (nPoints, 2 * nPoints))
|
|
|
|
# write indices
|
|
print("Writing indices to VTK file...")
|
|
for index in range(nPoints):
|
|
outfile.write('%10d %10d\n' % (1, index))
|
|
|
|
outfile.write('POINT_DATA %15d\n' % (nPoints))
|
|
outfile.write('SCALARS shotnumbers int %d\n' % (1))
|
|
outfile.write('LOOKUP_TABLE default\n')
|
|
|
|
# write shotnumber
|
|
print("Writing shotnumbers to VTK file...")
|
|
for shotnumber in shotnumbers:
|
|
outfile.write('%10d\n' % shotnumber)
|
|
|
|
outfile.close()
|
|
print("Wrote %d sources 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)
|