jpjones76 / SeisIO.jl

Julia language support for geophysical time series data
http://seisio.readthedocs.org
Other
47 stars 21 forks source link

merge! between SeisData, SeisChannel and SeisEvent, EventTraceData, EventChannel #59

Closed tclements closed 3 years ago

tclements commented 4 years ago

Starting a discussion on merging SeisData and SeisChanel to SeisEvent, EventTraceData and EventChannel in the Quake module. For the dataset with which I'm working, EQMag, EQLoc and SeisPha are stored separately from the seismic traces.

The workflow I'm currently using for merging SeisEvent with phase, location, and magnitude information and trace data

  1. Populate an EQLoc
  2. Populate an EQMag
  3. Create an empty EventTraceData
  4. Fill EventTraceData with EventChannels
  5. Populate a SeisHdr with EQLoc and EQMag
  6. Populate a SeisSrc
  7. Populate SeisEvent with SeisHdr, SeisSrc and EventTraceData
  8. Read waveform data into SeisData
  9. Merge .x, .t, .fs, .gain channels in SeisChannel into SeisEvent.data.EventChannel

Conversion from SeisData to EventTraceData currently works

julia> EventTraceData(SeisData())
EventTraceData with 0 channels (0 shown)
...

but conversion from SeisChannel to EventChannel is not yet implemented.

julia> EventChannel(SeisChannel())
ERROR: MethodError: no method matching EventChannel(::SeisChannel)

It would be nice to have a merge! function for SeisChannel and EventChannel which checks if each data field for common channels is empty. e.g. the SeisChannel will have non-empty .fs,.t and .x fields and the EventChannel will have non-empty .loc, .baz, .dist and possible .loc.

merge! with SeisData and EventTraceData could work similarly.

I'm working on a solution right now but happy to discuss if these additions are helpful/needed.

jpjones76 commented 4 years ago

I really like this idea. If you don't have time to crunch through a solution, I can start on this later in the week.

tclements commented 4 years ago

This is what I'm working with right now. The logic is to favor any (non-empty) data in SeisChannel and SeisData over anything in the SeisEvent. I could add a merge! method with SeisEvent and EventTraceData and SeisChannel /SeisData if we add matching + or append methods.

using SeisIO
using SeisIO.Quake
import SeisIO: merge, merge!

function merge!(EC::EventChannel,SC::SeisChannel)
    if EC.id != SC.id 
        throw(ErrorException("EventChannel ID $(EC.id) and SeisChannel ID $(SC.id) do not match."))
    end
    for f in SeisIO.datafields
        fval = getfield(SC,f)
        if !isempty(fval)
            setfield!(EC,f,fval)
        end
    end
    return nothing
end
merge!(SC::SeisChannel,EC::EventChannel,) = merge!(EC,SC)
merge(EC::EventChannel,SC::SeisChannel) = (U = deepcopy(EC);merge!(U,SC);return U)
merge(SC::SeisChannel,EC::EventChannel) = merge(EC,SC)

# this method should call merge!(SE,SD) once implemented
function merge(SE::SeisEvent,SD::SeisData)
    newEvent = SeisEvent()
    newEvent.hdr = SE.hdr
    newEvent.source = SE.source
    newEvent.data = merge(SE.data,SD)
    return newEvent
end
merge(SD::SeisData,SE::SeisEvent) = merge(SE,SD)

# this method should call merge!(ETD,SD) once implemented
# this currently only merges events with seismic data 
# need to check that SD has data 
function merge(ETD::EventTraceData,SD::SeisData)
    N = length(SD.id)
    ETDdict = Dict(ETD.id[ii]=>ii for ii = 1:length(ETD.id))
    ETDkeys = keys(ETDdict)
    newETD = EventTraceData(N)
    for ii = 1:N
        if in(SD.id[ii],ETDkeys)
            ind = ETDdict[SD.id[ii]]
            newETD[ii] = merge(SD[ii],ETD[ind])
        else
            newETD[ii] = EventTraceData(SD[ii])[1] # need to implement EventChannel(SeisChannel())
        end
    end
    return newETD 
end
merge(SD::SeisData,ETD::EventTraceData) = merge(ETD,SD)
jpjones76 commented 3 years ago

I think this can be implemented from what you've done here. Sorry for dropping the ball on this last month; several wildfires happened, followed by a long-distance move that couldn't be delayed. I'm back at it now.

jpjones76 commented 3 years ago

I'm still working on this. Fundamentally, for any of these methods to be plausible, there must be a way to merge! into a (single-channel) GphysChannel subtype, without converting to/from the multichannel equivalent. That's doable, but not straightforward.

At issue is that non-data fields must be merged in a way that's consistent with merge! for SeisData. That operation retains field values from the segment with the latest start time. If I want strict consistency, then some merge! operations need to move to helper functions. I'm working on that now.

jpjones76 commented 3 years ago

I need to ask a very basic question here, which I can't guess from context.

Do you actually merge your :x and :t fields at all, or are you merely copying them from a populated SeisChannel to the corresponding field of an empty EventChannel?

tclements commented 3 years ago

I'm just copying rather than merging :x or :t fields for now. I think we should do all mergeing on the SeisChannel <-> EventChannel level though.

jpjones76 commented 3 years ago

I've added GphysChannel merges in SeisIO v1.2.0. They use a partial match function that you could import and implement in your own code. If that matches what you need, please let me know, and I can close this issue.

tclements commented 3 years ago

Sounds good! I can close.