vpick script and rdsegy.m Matlab function
#! /bin/csh -f
setenv DISPLAY walrus:0.0 # CHANGE THIS FOR YOUR MACHINE!
if( $#argv < 2 ) then
echo "Usage: vpick filename start_rpnum [ rpnum_inc end_rpnum]"
exit 1
endif
set FILE = $1
set MATFILE = vpick.mat # must terminate with .mat for Matlab
set VPICKFILE = $FILE.vpick
set START_RPNUM = $2
if( $#argv < 3 ) then
set RPNUM_INC = 1
else
set RPNUM_INC = $3
endif
if( $#argv < 4 ) then
set END_RPNUM = $START_RPNUM
else
set END_RPNUM = $4
endif
set rpnum = $START_RPNUM
while ( $rpnum <= $END_RPNUM )
sioseis << eof
noecho procs diskin filter agc plot velan prout END
diskin
secs 4 fno $rpnum lno $rpnum
ipath $FILE END END
header # chnage the rp number into the water depth
fno 0 lno 99999 ftr 0 ltr 9999
l6 = r54
end
end
filter pass 20 220 dbdrop 48 minpha no end end
agc
winlen .5 END END
prout
fno 0 lno 99999 ftr 0 ltr 99999 END
END
velan
vels 1400 10 1800 nrp 1 type spec winlen .048 opath $MATFILE
END
END
plot
nibs 75 stime 3. nsecs 2.5 scalar -1 trpin 10 vscale 5. def 0.06
srpath sunfil ftag 1 taginc 1000 ann fanno fanno $rpnum END
END
END
eof
xloadimage -r 90 sunfil &
# Use Matlab Version 4.2c (Nov 23 1994) or newer or else change
# getline to ginput, which does not have a line connecting the picks.
matlab << eof1
load $MATFILE;
n = size(vel);nt = n(1);
nv = n(2);
rpno = vel(1,1);
st = vel(2,1);
dt = vel(3,1);
sv = vel(4,1);
dv = vel(5,1);
vel(1,1) = 0;vel(2,1) = 0;vel(3,1) = 0;vel(4,1) = 0;vel(5,1) = 0;
x = sv:dv:sv+nv*dv-dv;
y = -st:-dt:-(st+nt*dt-dt);
contour(x,y,vel,20)
title 'rp $rpnum';
xx=[];
yy=[];
n = 0;
but = 1;
while but == 1
[xi,yi,button] = ginput(1);
n = n + 1;
xx(n,1) = xi;
yy(n,1) = yi;
if button > 1
but = 0;
end
end
pause
n = size(xx);
fprintf('$VPICKFILE',' fno %.0f vtp ',$rpnum)
for i = 1:n(1)
fprintf('$VPICKFILE','%.0f %.3f ',xx(i),-yy(i))
end
fprintf('$VPICKFILE','end\n')
quit
eof1
@ rpnum = $rpnum + $RPNUM_INC
end
rm $MATFILE
function [seis,ntr,fmt,dsrt,dom,x,t0,nsamps,dt,nx,sh, ...
tr,rp,trp,trcid]=readsegy(filenm,a1,a2);
%
% Reads an SEG-Y trace into "ans"
% arguments:
% filename - The name of the SEG-Y file to be read
% a1 - The shot/rp number of the trace to be read
% a2 - The trace number to be read.
%
% Short program-driver for reading SEGY files
% The SEGY binary tape header parameters passed here are
% (all of them of 1x1 size, because it is a tape header,
% so you have to write it out just once):
% seis(nsamps by nx) Seismogram data
% ntr : Number of traces per shot (or record)
% fmt : SEGY format type, >4, "host" floating point
% dsrt : How is data sorted. 0 or 1 by shots; 2 by CDP gathers
% dom : The domain of the data (0 or 1 - time, 7 - tau-p, etc.)
% Variables are in the two temporary arrays -
% lbuf holds long integers and ibuf holds short integers
% The header parameters passed here are:
% x (nx by 1) : Range in km
% t0(nx by 1) : Start time in sec
%%% id(nx by 1) : Trace index integer (ie 1 through nx)
% nsamps(nx by 1) : Number of samples per trace
% dt(1 by 1) : Sample interval in sec/sample
% nx : Number of traces written in
% sh : Shot number
% tr : Trace within the shot
% rp : RP or CDP number
% trp : Trace number within RP or CDP
% trcfid : Trace ID; live = 1
fno = str2num(a1);
ftr = str2num(a2);
% Open file and read the "tape" header
fid = fopen(filenm, 'r');
if fid == -1
error('error opening SEG-Y file')
end
% Read the EBCDIC header, putting into temp arrays
% Useful information for subsequent writing is stored in appropriate variables
thbuf1 = fread(fid,3200)';
%disp( 'EBCDIC tape header ')
thbuf1(1,1);
%fseek(fid,3200,0);
thbuf2 = fread(fid,200,'short')';
ntr = thbuf2(7);
fmt = thbuf2(13);
dsrt = thbuf2(15);
dom = thbuf2(31);
tst = thbuf2(200);
if (fmt == 1)
exit('Cannot read IBM floating point format');
end
% Next read the trace headers and data, putting into temp arrays
% Reads each header twice, first putting into ibuf in short int format
% then into lbuf in long int format
% Array ibuf is (nx by 120) and array lbuf is (nx by 60)
buftr = 40; % Size of buffer increments
nx = 0; % Number of traces read
nread = 0;
doit = 1;
while (doit)
% T1 = fread(fid,120,'short')';
T1 = fread(fid,120,'int16')';
if (feof(fid) == 1) break; end
nread = nread + 1;
if (nread == 1 ) % Preassign Buffer space for trace headers/data
nsamps = T1(58);
bibuf = zeros(buftr,120);
blbuf = zeros(buftr, 60);
bseis = zeros(buftr,nsamps);
end
nsamps = T1(58);
% fseek(fid,-240,0)
fseek(fid,-240,'cof');
% T2 = fread(fid,60,'long')';
T2 = fread(fid,60,'int32')';
if (feof(fid) == 1)
fprintf(1,'reread of header failed.\n');
break;
end
if( fmt == 2 ) T3 = fread(fid,nsamps,'long')'; end;
if( fmt == 3 ) T3 = fread(fid,nsamps,'short')'; end;
if( fmt > 4 ) T3 = fread(fid,nsamps,'float')'; end;
if (feof(fid) == 1)
fprintf(1,'read of data failed.\n');
break;
end
nx = nx + 1;
bibuf(nx,:) = T1;
blbuf(nx,:) = T2;
if( blbuf(nx,7) == 0 )
no = blbuf(nx,3);
tr = blbuf(nx,4);
else
no = blbuf(nx,6);
tr = blbuf(nx,7);
end
bseis(nx,:) = T3;
if( no == fno )
if( tr == ftr )
doit = 0;
else
nx = nx - 1;
end
else
nx = nx - 1;
end
% Allocate additional storage space for data arrays
if ( nx > 0 )
if( rem(nx,buftr) == 0)
seis = [seis;bseis];
ibuf = [ibuf;bibuf];
lbuf = [lbuf;blbuf];
end
end
end
if ( nx == 0 )
quit
end
seis = bseis(1:nx,:);
ibuf = bibuf(1:nx,:);
lbuf = blbuf(1:nx,:);
fclose(fid);
% identify the variables inside the ibuf - the short integers
% or lbuf - the long integers, arrays to pass on
dt = bibuf(:,59) * 0.000001;
nsamps = bibuf(:,58);
x = blbuf(:,10) * 0.001;
t0 = bibuf(:,55) * 0.001;
sh = blbuf(:,3);
tr = blbuf(:,4);
rp = blbuf(:,6);
trp = blbuf(:,7);
trcid = bibuf(:,15);
%id = [:]'
seis = bseis(:,:)';
Return to EW0008 example.
Return to SIOSEIS examples.
Go to the list of seismic processes. Go to SIOSEIS introduction.