-
Notifications
You must be signed in to change notification settings - Fork 50
/
SeisWindow.jl
116 lines (102 loc) · 3.15 KB
/
SeisWindow.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
"""
SeisWindow(in,out;<keyword arguments>)
Window a seis file using header words.
# Arguments
* `in::AbstractString`: filename of input
* `out::AbstractString`: filename of output
# Keyword arguments
* `key`
* `minval`
* `maxval`
note that windowing along the time axis is achieved by using the key "t".
*Credits: AS, 2015*
"""
function SeisWindow(in::AbstractString,out::AbstractString;key=[],minval=[],maxval=[])
DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
extent = ReadTextHeader(in)
tmin = extent.o1
tmax = extent.o1 + extent.d1*(extent.n1-1)
for ikey=1:length(key)
if key[ikey] == "t"
tmin = minval[ikey]
tmax = maxval[ikey]
println("tmin= ",tmin," tmax= ",tmax)
end
end
#itmin is the sample at which tmin is saved
itmin = convert(Int32,round((tmin - extent.o1)/extent.d1) + 1)
if (itmin < 1)
itmin = 1
end
itmax = convert(Int32,round((tmax - extent.o1)/extent.d1))
if (itmax) > extent.n1
itmax = extent.n1
end
println("itmin= ",itmin," itmax= ",itmax)
SeisWindowHeaders(in, out; key=key, minval=minval, maxval=maxval, tmin=tmin,
tmax=tmax)
FetchTraces(in,out;itmin=itmin,itmax=itmax)
tmp = join(["tmp_SeisWindow_",string(round(Int,rand()*100000))])
@compat SeisProcessHeaders(out, tmp, [UpdateHeader],
[Dict(:itmin=>itmin,:itmax=>itmax)])
filename_h_tmp = join([DATAPATH tmp "@headers@"])
filename_h_out = join([DATAPATH out "@headers@"])
cp(filename_h_tmp,filename_h_out,remove_destination=true);
rm(filename_h_tmp);
end
function FetchTraces(in::AbstractString, out::AbstractString; ntrace=500, itmin=round(Int,1),
itmax=round(Int,9999999999))
NX = GetNumTraces(out)
println("NX in FetchTraces= ",NX)
itrace = round(Int,1)
filename_data_in = ParseDataName(in)
DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
filename_data_out = join([DATAPATH out "@data@"])
stream_in = open(filename_data_in)
stream_out = open(filename_data_out,"w")
extent = ReadTextHeader(in)
nt = extent.n1
if itmin < 1
itmin = round(Int,1)
end
if itmax > nt
itmax = round(Int,nt)
end
while itrace <= NX
if (itrace > 1)
stream_out = open(filename_data_out,"a")
end
nx = NX - itrace + 1
ntrace = nx > ntrace ? ntrace : nx
h = SeisReadHeaders(out,group="some",itrace=itrace,ntrace=ntrace)
d = zeros(Float32,itmax-itmin+1,ntrace)
SeekTraces!(d,stream_in,h,itmin,itmax,nt,round(Int,ntrace))
write(stream_out,d)
close(stream_out)
itrace += ntrace
end
close(stream_in)
end
function SeekTraces!{T}(d::AbstractArray{T,2}, stream_in::IOStream,
h::Array{Header,1},itmin,itmax,nt,ntrace)
d1 = zeros(Float32,nt)
for ix = 1 : ntrace
position = 4*nt*(h[ix].tracenum-1)
seek(stream_in,position)
d1 = read(stream_in,Float32,nt)
for it = itmin : itmax
d[it - itmin + 1,ix] = d1[it]
end
end
nothing
end
function UpdateHeader(h;itrace=1,itmin=1,itmax=9e9)
ot = (itmin-1)*h[1].d1 + h[1].o1
nt = itmax - itmin + 1
for j = 1 : length(h)
h[j].o1 = ot
h[j].n1 = nt
h[j].tracenum = j + itrace - 1
end
return h
end