Dev (#84)

* update smoothing with edge smoothing (#55)

* Add NodalGPU to dev  (#56)

* Update Nodal stuff

* Fix CorrData filtering typo

* Fix wrong id when file starts before midnight

* Typo in RawData fix

* Fix shorten bug with newlag/sampling rate mismatch

* Allow whiten! to work with ComplexF64

* Remove Plots and DataFrames dependencies!

Project.toml

@@ -7,30 +7,26 @@ version = "0.5.0"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SeisIO = "b372bb87-02dd-52bb-bcf6-c30dd83fd342"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Adapt = "1.0.1, 2.0, 3.0"
-CUDA = "1.0"
 DSP = "0.6"
-DataFrames = "0.20, 0.21"
 FFTW = "1.1"
-GPUArrays = "3.1, 4.0, 5.0"
+GPUArrays = "3.1, 4.0, 5.0, 6.2"
 Glob = "1.2"
 JLD2 = "0.1, 0.2"
-Plots = "0.29, 1.0"
 SeisIO = "1.0"
 StatsBase = "0.32,0.33"
 julia = "1"

src/ArrayFuncs.jl

@@ -33,6 +33,8 @@ detrend!(R::RawData) = detrend!(R.x)
 detrend(R::RawData) = (U = deepcopy(R); detrend!(U.x); return U)
 detrend!(C::CorrData) = detrend!(C.corr)
 detrend(C::CorrData) = (U = deepcopy(C); detrend!(U.corr); return U)
+detrend!(N::NodalData) = detrend!(N.data)
+detrend(N::NodalData) = (U = deepcopy(N); detrend!(U.data); return U)
 
 
 """
@@ -51,6 +53,8 @@ demean!(R::RawData) = demean!(R.x)
 demean(R::RawData) = (U = deepcopy(R); demean!(U.x); return U)
 demean!(C::CorrData) = demean!(C.corr)
 demean(C::CorrData) = (U = deepcopy(C); demean!(U.corr); return U)
+demean!(N::NodalData) = demean!(N.data)
+demean(N::NodalData) = (U = deepcopy(N); demean!(U.data); return U)
 
 """
  taper!(A,fs; max_percentage=0.05, max_length=20.)
@@ -92,6 +96,12 @@ taper!(C::CorrData; max_percentage::AbstractFloat=0.05,
 taper(C::CorrData; max_percentage::AbstractFloat=0.05,
  max_length::Real=20.) = (U = deepcopy(C); taper!(U.corr,U.fs,
  max_percentage=max_percentage,max_length=max_length); return U)
+taper!(N::NodalData; max_percentage::AbstractFloat=0.05,
+ max_length::Real=20.) = taper!(N.data,N.fs[1],max_percentage=max_percentage,
+ max_length=max_length)
+taper(N::NodalData; max_percentage::AbstractFloat=0.05,
+ max_length::Real=20.) = (U = deepcopy(N); taper!(U.data,U.fs[1],
+ max_percentage=max_percentage,max_length=max_length); return U)
 
 
 """

src/Plotting/plotting.jl

@@ -1,17 +1,33 @@
-export corrplot
+@recipe function f(C::CorrData)
+ # set up the subplots
+ grid := false
+ layout := (2,1)
 
-"""
-
- corrplot(C)
-
-Plot a CorrData matrix.
-"""
-function corrplot(C::CorrData)
  lags = -C.maxlag:1/C.fs:C.maxlag
  times = Dates.format.(Dates.unix2datetime.(C.t),"yyyy/m/d HH:MM")
- Cstack = stack(C,allstack=true)
- plot(
- heatmap(lags,times,C.corr',c=:balance,legend=:none),
- plot(lags,Cstack.corr,c=:black,linewidth=1.5,legend=:none,xlabel="Lag [s]"),
- layout = grid(2,1,heights=[0.75,0.25]), link=:x,dpi=1000)
-end
+
+ # main heatmap
+ @series begin
+ legend := false
+ seriestype := :heatmap
+ ytickfontsize --> 10 
+ xtickfontsize --> 8
+ subplot := 1
+ seriescolor --> :balance
+ lags,times,C.corr'
+ end
+
+ # bottom stack
+ @series begin
+ subplot := 2
+ label := false
+ title --> C.name
+ xguide --> "Lag [s]"
+ titlefontsize --> 6 
+ ytickfontsize --> 10 
+ xtickfontsize --> 8
+ seriescolor --> :black
+ xlims --> (lags[1],lags[end])
+ lags,stack(C,allstack=true).corr
+ end
+end

src/SeisNoise.jl

@@ -1,8 +1,7 @@
 module SeisNoise
 
-using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO
-using Statistics, StatsBase, Plots, Distributed
-using CUDA, Adapt, GPUArrays
+using Dates, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO, SeisIO.Nodal
+using Statistics, StatsBase, Distributed, CUDA, Adapt, GPUArrays, RecipesBase
 
 # check use of cuda
 const use_cuda = Ref(false)
@@ -13,6 +12,7 @@ end
 
 # import types first
 include("Types/NoiseData.jl")
+include("Types/NodalData.jl")
 include("Types/show.jl")
 
 # import pre and post processing tools
@@ -28,6 +28,7 @@ include("stacking.jl")
 # import routines for doin' stuff
 include("compute_fft.jl")
 include("correlation.jl")
+include("nodalcorrelation.jl")
 include("rotation.jl")
 include("Plotting/plotting.jl")
 

src/Types/NodalData.jl

@@ -0,0 +1,80 @@
+export gpu,cpu, NodalCorr
+
+mutable struct NodalCorr <: NoiseData
+ n::Int64
+ id::Array{String,1} # id
+ name::Array{String,1} # name
+ loc::Array{<:InstrumentPosition,1} # loc
+ fs::Float64 # fs
+ src::String # src
+ t::Array{Array{Int64,2},1} # time
+ corr::AbstractArray{Float32, 2} # actual data
+
+ function NodalCorr(
+ n ::Int64,
+ id ::Array{String,1},
+ name ::Array{String,1},
+ loc ::Array{<:InstrumentPosition,1},
+ fs ::Float64,
+ src ::String,
+ t ::Array{Array{Int64,2},1} ,
+ corr ::AbstractArray{<:AbstractFloat,2},
+ )
+
+ return new(n, id, name, loc, fs, src, t, corr)
+ end
+end
+
+NodalCorr(;
+ n ::Int64 = 0,
+ id ::Array{String,1} = Array{String}(undef,0),
+ name ::Array{String,1} = Array{String}(undef,0),
+ loc ::Array{InstrumentPosition,1} = Array{InstrumentPosition}(undef,0), 
+ fs ::Float64 = zero(Float64),
+ src ::String = "",
+ t ::Array{Array{Int64,2},1} = Array{Array{Int64,2},1}(undef,0),
+ corr ::AbstractArray{<:AbstractFloat,2} = Array{Float32,2}(undef, 0, 2),
+ ) = CorrData(n, id, name, loc, fs, src, t, corr)
+
+
+# functions for adapting NodalData to GPU 
+cpu(m::NodalData) = adapt(Array,m)
+gpu(x::NodalData) = use_cuda[] ? cu(x) : x
+
+# function for adapting NodalCorr to GPU 
+cpu(m::NodalCorr) = adapt(Array,m)
+gpu(x::NodalCorr) = use_cuda[] ? cu(x) : x
+
+Adapt.adapt_structure(to, N::NodalData) = NodalData(
+ N.n,
+ N.ox,
+ N.oy,
+ N.oz,
+ N.info,
+ N.id,
+ N.name,
+ N.loc,
+ N.fs,
+ N.gain,
+ N.resp,
+ N.units,
+ N.src,
+ N.misc,
+ N.notes,
+ N.t,
+ adapt(to,N.data),
+)
+
+Adapt.adapt_structure(to, N::NodalCorr) = NodalCorr(
+ N.n,
+ N.id,
+ N.name,
+ N.loc,
+ N.fs,
+ N.gain,
+ N.resp,
+ N.units,
+ N.src,
+ N.t,
+ adapt(to,N.corr),
+)

src/Types/NoiseData.jl

@@ -105,7 +105,7 @@ mutable struct RawData <: NoiseData
  # get starting and ending indices
  startind, endind = slide_ind(startslice, endslice, S.fs[1], S.t[1])
  x, starts = slide(@view(S.x[1][startind:endind]), cc_len, cc_step, S.fs[1], startslice,endslice)
- return new(S[1].id,Dates.format(u2d(S.t[1][1,2]*1e-6),"Y-mm-dd"),S.loc[1],
+ return new(S[1].id,Dates.format(u2d(starts[1]),"Y-mm-dd"),S.loc[1],
  S.fs[1],S.gain[1],1. / cc_len,S[1].fs/2,cc_len,cc_step,false,
  "",S[1].resp,S.misc[1],S.notes[1],starts,x)
  end
@@ -140,7 +140,7 @@ mutable struct RawData <: NoiseData
  # get starting and ending indices
  startind, endind = slide_ind(startslice, endslice, C.fs,C.t)
  x, starts = slide(@view(C.x[startind:endind]), cc_len, cc_step, C.fs, startslice,endslice)
- return new(C.id,Dates.format(u2d(C.t[1,2]*1e-6),"Y-mm-dd"),C.loc,C.fs,
+ return new(C.id,Dates.format(u2d(starts[1]),"Y-mm-dd"),C.loc,C.fs,
  C.gain,1. / cc_len,C.fs/2,cc_len,cc_step,false,"", C.resp,
  C.misc,C.notes,starts,x)
  end

src/correlation.jl

@@ -152,22 +152,21 @@ function correlate(FFT1::FFTData, FFT2::FFTData, maxlag::Real;corr_type::String=
 end
 
 """
- whiten!(A, freqmin, freqmax, fs, pad=50)
+ whiten!(A, freqmin, freqmax, fs, N; pad=50)
 
 Whiten spectrum of rfft `A` between frequencies `freqmin` and `freqmax`.
 Returns the whitened rfft of the time series.
 
 # Arguments
 - `A::AbstractArray`: Time series.
-- `fs::Real`: Sampling rate of time series `A`.
 - `freqmin::Real`: Pass band low corner frequency.
 - `freqmax::Real`: Pass band high corner frequency.
+- `fs::Real`: Sampling rate of time series `A`.
 - `N::Int`: Number of input time domain samples for each rfft.
 - `pad::Int`: Number of tapering points outside whitening band.
 """
-function whiten!(A::AbstractArray{Complex{Float32}}, freqmin::Real,
- freqmax::Real, fs::Real,N::Int;pad::Int=50)
- T = real(eltype(A))
+function whiten!(A::AbstractArray{Complex{T}}, freqmin::Real,
+ freqmax::Real, fs::Real,N::Int;pad::Int=50) where T <: AbstractFloat
  Nrows,Ncols = size(A)
 
  # get whitening frequencies
@@ -268,6 +267,18 @@ end
 whiten(R::RawData,freqmin::Real, freqmax::Real; pad::Int=50) =
  (U = deepcopy(R); whiten!(U,freqmin,freqmax,pad=pad);return U)
 
+function whiten!(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50)
+ @assert freqmin > 0 "Whitening frequency must be greater than zero."
+ @assert freqmax <= N.fs[1] / 2 "Whitening frequency must be less than or equal to Nyquist frequency."
+ Npts = size(N.data,1) # number of data points
+ FFT = rfft(N.data,1)
+ whiten!(FFT,freqmin,freqmax,N.fs[1], Npts, pad=pad)
+ N.data .= irfft(FFT,Npts,1)
+ return nothing
+end
+whiten(N::NodalData,freqmin::Real, freqmax::Real; pad::Int=50) =
+ (U = deepcopy(N); whiten!(U,freqmin,freqmax,pad=pad);return U)
+
 """
 
  coherence!(F,half_win, water_level)

src/deprecated.jl

@@ -1,8 +1,12 @@
-export compute_fft, compute_cc
+export compute_fft, compute_cc, corrplot
 
 @deprecate compute_fft() rfft()
 @deprecate compute_cc() correlate()
 
+function corrplot(C::CorrData)
+ @warn "corrplot has been deprecated. Please use: `using Plots;plot(C)` to plot CorrData."
+end
+
 """
 
  compute_fft(R)

src/distance.jl

@@ -374,3 +374,22 @@ function get_loc(C::CorrData)
 
  return C.loc, loc2
 end
+
+function get_azi(loc1::NodalLoc,loc2::NodalLoc)
+ azi = atan(loc2.y - loc1.y, loc2.x - loc1.x) * 180 / π 
+ if azi <= 90 
+ azi = 90 - azi 
+ else
+ azi = 450 - azi 
+ end
+ return azi 
+end
+
+function get_azi(N::NodalData)
+ out = zeros(N.n)
+ for ii = 1:N.n-1
+ out[ii] = get_azi(N.loc[ii],N.loc[ii+1])
+ end
+ out[end] = out[end-1]
+ return out
+end

src/filter.jl

@@ -72,9 +72,9 @@ bandpass(R::RawData,freqmin::Real,freqmax::Real;
  setfield!(U,:freqmin,Float64(freqmin));
  setfield!(U,:freqmax,Float64(min(freqmax,U.fs/2)));return U)
 bandpass!(C::CorrData,freqmin::Real,freqmax::Real;
-  corners::Int=4,zerophase::Bool=true) = (bandpass!(C.corr,freqmin,freqmax,
-  C.fs,corners=corners,zerophase=zerophase);setfield!(C,:freqmin,Float64(freqmin));
-  setfield!(C,:freqmax,Float64(min(freqmax,C.fs/2)));return nothing)
+ corners::Int=4,zerophase::Bool=true) = (bandpass!(C.corr,freqmin,freqmax,
+ C.fs,corners=corners,zerophase=zerophase);setfield!(C,:freqmin,Float64(freqmin));
+ setfield!(C,:freqmax,Float64(min(freqmax,C.fs/2)));return nothing)
 bandpass(C::CorrData,freqmin::Real,freqmax::Real;
  corners::Int=4,zerophase::Bool=true) = (U = deepcopy(C);bandpass!(U.corr,
  freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);
@@ -91,6 +91,13 @@ bandpass!(S::SeisData, freqmin::Real, freqmax::Real;
  fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
 bandpass(S::SeisData, freqmin::Real, freqmax::Real; corners::Int=4,
  zerophase::Bool=true) = filtfilt(S,fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandpass")
+bandpass!(N::NodalData,freqmin::Real,freqmax::Real;
+ corners::Int=4,zerophase::Bool=true) = (bandpass!(N.data,freqmin,freqmax,
+ N.fs[1],corners=corners,zerophase=zerophase);return nothing)
+bandpass(N::NodalData,freqmin::Real,freqmax::Real;
+ corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandpass!(N.data,
+ freqmin,freqmax,N.fs[1],corners=corners,zerophase=zerophase);
+ return U)
 
 """
  bandstop!(A,freqmin,freqmax,fs,corners=4,zerophase=true)
@@ -175,6 +182,12 @@ bandstop!(S::SeisData, freqmin::Real, freqmax::Real;
 bandstop(S::SeisData,freqmin::Real, freqmax::Real; corners::Int=4,
  zerophase::Bool=true) = filtfilt(S,
  fl=Float64(freqmin),fh=Float64(freqmax),np=corners,rt="Bandstop")
+bandstop!(N::NodalData,freqmin::Real,freqmax::Real;
+ corners::Int=4,zerophase::Bool=true) = bandstop!(N.data,freqmin,freqmax,
+ N.fs[1],corners=corners,zerophase=zerophase)
+bandstop(N::NodalData,freqmin::Real,freqmax::Real;
+ corners::Int=4,zerophase::Bool=true) = (U = deepcopy(N);bandstop!(N.data,
+ freqmin,freqmax,U.fs[1],corners=corners,zerophase=zerophase);return U)
 
 """
 lowpass(A,freq,fs,corners=4,zerophase=true)
@@ -250,6 +263,12 @@ lowpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
  filtfilt!(S,fh=Float64(freq),np=corners,rt="Lowpass")
 lowpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
  filtfilt(S,fh=Float64(freq),np=corners,rt="Lowpass")
+lowpass!(N::NodalData,freq::Real; corners::Int=4,
+ zerophase::Bool=true) = (lowpass!(N.data,freq,N.fs[1],corners=corners,
+ zerophase=zerophase);return nothing)
+lowpass(N::NodalData,freq::Real; corners::Int=4,
+ zerophase::Bool=true) = (U = deepcopy(N);lowpass!(N.data,freq,U.fs[1],
+ corners=corners,zerophase=zerophase);return U)
 
 """
  highpass(A,freq,fs,corners=4,zerophase=true)
@@ -321,6 +340,12 @@ highpass!(S::SeisData, freq::Real;corners::Int=4, zerophase::Bool=true) =
  filtfilt!(S,fl=Float64(freq),np=corners,rt="Highpass")
 highpass(S::SeisData,freq::Real; corners::Int=4,zerophase::Bool=true) =
  filtfilt(S,fl=Float64(freq),np=corners,rt="Highpass")
+highpass!(N::NodalData,freq::Real; corners::Int=4,
+ zerophase::Bool=true) = (highpass!(N.data,freq,N.fs[1],corners=corners,
+ zerophase=zerophase);return nothing)
+highpass(N::NodalData,freq::Real; corners::Int=4,
+ zerophase::Bool=true) = (U = deepcopy(N);highpass!(N.data,freq,U.fs[1],
+ corners=corners,zerophase=zerophase);return U)
 
 """
  envelope(A)