Issues getting synthetics from FK database

[ammcpherson]

Good afternoon,

I am trying to save synthetics from MTUQ to visualize against real data. I have done this with SPECFEM-generated Green's functions, but when I attempt to do the same thing using an FK database, no synthetics are generated. I cannot share the FK database of course, but these are the basic commands of what I am doing:

fk_path = '/store/wf/FK_synthetics/'+model
db = open_db(fk_path,format='FK')

origin = Origin({
        'time': '2022-09-24T15:18:54.694000Z',
        'latitude': 61.4915,
        'longitude': -145.5887,
        'depth_in_m': 41000.0,
        'id': str(eid)
})

data = read(path_data, format='sac', 
            event_id=event_id,
            station_id_list=station_id_list,
            tags=['units:m', 'type:velocity']) 

data.sort_by_azimuth()
stations = data.get_stations()

greens = db.get_greens_tensors(stations, origin)

synthetics = greens.get_synthetics(MomentTensor([269000000000000,-15643000000000000,15374000000000000,5999000000000000,7576000000000000,4205000000000000]),components=['Z','R','T'])

os.makedirs(eid+'_synthetics',exist_ok=True)
synthetics.write(eid+'_synthetics/%s/' % model,format='sac') 

As I said, this works when I am using SPECFEM-generated GFs, and the only thing that changes is how I open the GF database:

SF_path = './GFs/%s/%s' % (eid,model)
db = open_db(SF_path,format='SPECFEM3D')

Is there something that needs to change about my approach to get synthetics from an FK database?

Thanks, Amanda

[rmodrak]

Hi Amanda, Thanks for the report. Maybe add some print() or type() statements to try to get to the bottom of it? I will take a look as well.

On Wed, Oct 23, 2024 at 5:26 PM Amanda McPherson @.***> wrote:

Good afternoon,

I am trying to save synthetics from MTUQ to visualize against real data. I have done this with SPECFEM-generated Green's functions, but when I attempt to do the same thing using an FK database, no synthetics are generated. I cannot share the FK database of course, but these are the basic commands of what I am doing:

fk_path = '/store/wf/FK_synthetics/'+model db = open_db(fk_path,format='FK')

origin = Origin({ 'time': '2022-09-24T15:18:54.694000Z', 'latitude': 61.4915, 'longitude': -145.5887, 'depth_in_m': 41000.0, 'id': str(eid) })

data = read(path_data, format='sac', event_id=event_id, station_id_list=station_id_list, tags=['units:m', 'type:velocity'])

data.sort_by_azimuth() stations = data.get_stations()

greens = db.get_greens_tensors(stations, origin)

synthetics = greens.get_synthetics(MomentTensor([269000000000000,-15643000000000000,15374000000000000,5999000000000000,7576000000000000,4205000000000000]),components=['Z','R','T'])

os.makedirs(eid+'_synthetics',exist_ok=True) synthetics.write(eid+'_synthetics/%s/' % model,format='sac')

As I said, this works when I am using SPECFEM-generated GFs, and the only thing that changes is how I open the GF database:

SF_path = './GFs/%s/%s' % (eid,model) db = open_db(SF_path,format='SPECFEM3D')

Is there something that needs to change about my approach to get synthetics from an FK database?

Thanks, Amanda

— Reply to this email directly, view it on GitHub https://github.com/uafgeotools/mtuq/issues/277, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABCGSSRYDZE7QDXFJQ3SRL3Z5AWD5AVCNFSM6AAAAABQP3QQNOVHI2DSMVQWIX3LMV43ASLTON2WKOZSGYYTAMBTGUZTQMI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[rmodrak]

It looks like the network, station, and location codes are not getting copied to the synthetics. The result is the following output file names. ...R.sac ...T.sac ...Z.sac

[rmodrak]

One way to get the correct filenames is to manually pass a list of ObsPy stats objects via the stats keyword argument when calling get_synthetics().

The syntax would be something like

synthetics = greens.get_synthetics(mt, stats=stations, components=components, mode='map')

This gets complicated though, I would have to look closely in order to recall the exact correct syntax.

Of course, we want these station codes to get passed seamlessly... trying to figure out a larger fix

[ammcpherson]

Thanks for your quick reply, @rmodrak!

I am having a little luck with a small code block along the lines of:

stats = {}
count = 0
for sta in stations:
    stats[count] = {}
    stats[count]['network'] = sta['network']
    stats[count]['station'] = sta['station']
    stats[count]['location'] = 'S3'
    for ii in range(3):
                if ii == 0:
                    stats[count]['channel'] = 'Z'
                elif ii == 1:
                    count += 1
                    stats[count] = stats[count-1]
                    stats[count]['channel'] = 'R'
                elif ii == 2:
                    count += 1
                    stats[count] = stats[count-1]
                    stats[count]['channel'] = 'T'     
    count += 1

In which I now get 1 file written out (there should be closer to 200 files) with the naming convention 'AK.PAX.S3.T.sac'. The 'S3' comes from the naming convention written out from the SPECFEM GF generated synthetics. This only works with mode='apply' in get_synthetics, though. When I try to use mode='map' I seem to encounter issues in _allocate_stream, where the function is making some stats object but also attempting to use the one I have passed to get_synthetics.

[rmodrak]

Very helpful troubleshooting, thanks.

If you needed a quick workaround, maybe try modifying the write() method: https://github.com/uafgeotools/mtuq/blob/231cb43044ae336b20134d942135aeca50ef4776/mtuq/dataset.py#L293

It should be possible to pass the stations list as an argument to the write() method, then use that to build the filenames.

[ammcpherson]

I think part of the problem now is that the FK GFs do not have station names added to the traces when MTUQ puts them all together into a stream - it appears that I would have to do some source/receiver distance calculation to then match that up to the proper set of traces to get the synthetics.

[ammcpherson]

Is there any reason that adding

trace.stats.station = station.station
trace.stats.location = station.location
trace.stats.network = station.network

to mtuq/io/clients/FK_SAC.py roughly here would not solve this problem entirely? Would it somehow affect the synthetics negatively?

Thanks.

[rmodrak]

Thanks for these comments-- agreed, there is probably some approach along the lines you mentioned.

[rmodrak]

When get_synthetics() is called, there is an option to create new objects to hold the synthetics. (This can be helpful, say, if we want to have two different sets of synthetics at the same time for comparison. As a sort of convention, mtuq and other better known Python packages use the inplace argument in connection with this.)

So on occasion, metadata will need to get copied once again to these newly created Stream or Trace objects. There's probably some way to copy the metadata from the streams holding the original Green's functions, as you suggest. The approach below I think is not too dissimilar. https://github.com/uafgeotools/mtuq/pull/279/commits/76688989e1f33357de2cd16a01a3117dba04c3c9