Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

gamma to NNLoc small snippet #259

Open
10go-L opened this issue Jan 3, 2024 · 1 comment
Open

gamma to NNLoc small snippet #259

10go-L opened this issue Jan 3, 2024 · 1 comment
Labels
enhancement New feature or request

Comments

@10go-L
Copy link

10go-L commented Jan 3, 2024

Hi, just sharing a small code to get the results from gamma to NNLoc:

catalogs, assignments = association(pick_df, station_df, config, method=config["method"])
catalog = pd.DataFrame(catalogs)
assignments = pd.DataFrame(assignments, columns=["pick_index", "event_index", "gamma_score"])
print('\ntotal number of earthquakes: ', len(catalog))

wgs84 = CRS.from_epsg(4326)
utm_zone_18s = CRS.from_epsg(32718)  # UTM zone 18S
transformer = Transformer.from_crs(utm_zone_18s, wgs84, always_xy=True)
catalog["longitude"], catalog["latitude"] = zip(*catalog.apply(
    lambda x: transformer.transform(x["x(km)"] * 1000, x["y(km)"] * 1000), axis=1))

catalog["depth_km"] = catalog["z(km)"]
times = mdates.date2num(catalog['time'])
# catalog.to_csv("gamma_events.csv", index=False, float_format="%.3f", date_format='%Y-%m-%dT%H:%M:%S.%f')

pick_df = pick_df.join(assignments.set_index("pick_index")).fillna(-1).astype({'event_index': int})

# Initialize a list to store strings for NLL phase input file
phase_file_str = []

# Filter out the rows where event_index is -1 and group the DataFrame by 'event_index'
for _, group in pick_df[pick_df['event_index'] != -1].groupby('event_index'):
    group.reset_index(drop=True, inplace=True)  # Reset indices for each dataframe for looping through

    # Iterate through each pick in the group
    for idx, pick in group.iterrows():
        if idx == 0 and phase_file_str:
            phase_file_str.append(" ")  # Add a blank line after each event
        
        # Construct each line for NLL phase file
        this_line = f"{pick['id'].ljust(5)}"  # Station
        this_line += "? ? ? "  # Instrument, Component, Onset
        this_line += f"{pick['type'].upper().ljust(4)}"  # Phase type (P or S), converted to uppercase
        this_line += "? "  # First motion
        this_line += pick['timestamp'].strftime("%Y%m%d %H%M %S.%f")[:20]  # Timestamp
        this_line += " GAU "  # Error type = Gaussian

        # Conditional GAU error based on 'prob' value
        if pick['prob'] >= 0.85:
            this_line += "5.00e-02 "  # 0.05 sec err for prob >= 0.85
        elif pick['prob'] >= 0.7:
            this_line += "1.00e-01 "  # 0.10 sec err for prob >= 0.7
        elif pick['prob'] >= 0.55:
            this_line += "2.00e-01 "  # 0.20 sec err for prob >= 0.55
        else:
            this_line += "3.33e-01 "  # 0.333 sec err for prob >= 0.4

        this_line += "-1.00e+00 -1.00e+00 -1.00e+00"  # Coda duration, amplitude, period (not used)

        phase_file_str.append(this_line)  # Append line to list

# Write new phase file
with open('./NLL/picks/pickfile_s.hpf', 'w') as f:
    for item in phase_file_str:
        f.write("%s\n" % item)
@yetinam yetinam added the enhancement New feature or request label Jan 8, 2024
@yetinam
Copy link
Member

yetinam commented Jan 8, 2024

Hi @10go-L ,

thanks for the example. Would you mind creating a PR to add this as a small script in the contrib directory? I think the contrib directory would be a good place to have this code snippet.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

2 participants