CC of horizontal components (rotate cc to obtain RR,TT)

[kjdavidson]

I would like to generate cross-correlations of the horizontal (Radial and Transverse) components from 3 component ZNE data. The documentation suggests that best practice is to compute the full tensor and then rotate afterwards. I'm looking for some guidance on how to do this last step please!

I believe the relevant configuration setting is:
$ msnoise config get components_to_compute M components_to_compute: ZZ,ZE,ZN,EZ,EE,EN,NZ,NE,NN

This gives me 3 output directories when I run compute_cc - ZE, ZN, ZZ.

I would like to take the .mseed files from each of these directories and rotate them in ObsPy to produce ZRT output, using the ZE file as my E component of my stream, ZN as my N component, etc. The backazimuth that I would rotate to is that between the first and second station in my pair.

Is this the correct approach?

Thanks in advance for any help!

[ThomasLecocq]

Hey, sorry for the delay. This code will be part of the next release indeed.

In principle, you have to compute EE,NN,EN,NE for computing RR,TT.

The following code will read the REF stacks and compute the radial-radial and transverse-transverse components, using the station database for defining the azimuths:

import os, glob, logging
from obspy import read

from msnoise.api import *
db = connect()

filterid = 2

c = np.cos
s = np.sin

for comp in ["ZZ","EE","NN","EN","NE","RR","TT"]:
    dir = "STACKS/%02i/REF/%s"%(filterid, comp)
    if not os.path.isdir(dir):
        os.makedirs(dir)

for sta1,sta2 in get_station_pairs(db):
    station1 = "%s.%s" % (sta1.net, sta1.sta)
    station2 = "%s.%s" % (sta2.net, sta2.sta)
    print(station1, station2)
    comps = {}
    for comp in ["ZZ","EE","NN","EN","NE"]:
        tmp1 = read("STACKS/%02i/REF/%s/%s_%s.MSEED"%(filterid, comp, station1.replace(".","_"),station2.replace(".","_")), format="MSEED")[0]
        comps[comp] = tmp1.data

    X0, Y0, c0 = (sta1.X, sta1.Y, sta1.coordinates)
    X1, Y1, c1 = (sta2.X, sta2.Y, sta1.coordinates)

    if c0 == c1:
        coordinates = c0
    else:
        coordinates = 'MIX'

    az = azimuth(coordinates, X0, Y0, X1, Y1)
    dist = get_interstation_distance(sta1,sta2)
    print("Azimuth", az)
    print("Distance", dist)
    t = np.deg2rad(az)
    p = np.deg2rad(az-180)
    rotation_operator = np.array([[-c(t)*c(p), c(t)*s(p),-s(t)*s(p), s(t)*c(p)],
                                                    [-s(t)*s(p),-s(t)*c(p),-c(t)*c(p),-c(t)*s(p)],
                                                    [-c(t)*s(p),-c(t)*c(p), s(t)*c(p), s(t)*s(p)],
                                                    [-s(t)*c(p), s(t)*s(p), c(t)*s(p),-c(t)*c(p)]])
    [TT, RR, TR, RT] = np.dot(rotation_operator, np.array([comps["EE"], comps["EN"], comps["NN"], comps["NE"]]))
    
    tmp1.data = TT
    tmp1.write("STACKS/%02i/REF/TT/%s_%s.MSEED"%(filterid, station1.replace(".","_"),station2.replace(".","_")), format="MSEED")
    tmp1.data = RR
    tmp1.write("STACKS/%02i/REF/RR/%s_%s.MSEED"%(filterid, station1.replace(".","_"),station2.replace(".","_")), format="MSEED")
    del tmp1
    

[kjdavidson]

Fantastic, thank you very much for the comprehensive reply!

[ThomasLecocq]

If you could make some short tests of it (I did, but it was a long time ago):
For example: for a small dataset:

Part 1:
msnoise config set components_to_compute=ZZ,RR,TT
msnoise compute_cc_rot
msnoise stack -r
mv STACKS/01 STACKS/91

Part 2:
msnoise config set components_to_compute=ZZ,EE,NN,EN,NE
msnoise compute_cc
msnoise stack -r
then run the rotation code

Part 3:
Compare the REFs in folder 01 and 91

That would be great !

[kjdavidson]

Results against my test dataset, the REF stacks generated are not precisely identical between each method. Here is an example from one station pair:

My dataset was using UTM coordinates, which I needed to convert to DEG to work with the rotation code, I'm not sure if this could explain the difference between the two REF stacks? I will re-run the CCs using the same coordinate system in both, but that takes some time:

compute_cc_rot took 9 hours, 9 minutes, 47 seconds
compute_cc took 3 hours, 42 minutes, 32 seconds
rotation code ran in a few seconds

Dataset is one year of data, 7 stations (6942 files, ~50GB). I used sqlite for the database and configured a single filter as per the 'simple msnoise howto' in the docs.

So a big improvement on processing time using compute_cc and the new rotation code!

[ThomasLecocq]

Hi, thanks for the test! They are indeed different, but that makes sense and it looks like the main features are all well reproduced ! Good! -- Note normally the coordinate DEG or UTM shouldn't make a difference, I'll double check!