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

Velest #18

Closed
wants to merge 11 commits into from
Closed

Velest #18

wants to merge 11 commits into from

Conversation

HaoZhang2021
Copy link

add velest exporter

Copy link
Member

@miili miili left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice! The code shrunk incredibly!

Please see my review suggestions, contact me if you have any questions.

src/qseek/exporters/velest.py Outdated Show resolved Hide resolved
src/qseek/exporters/velest.py Outdated Show resolved Hide resolved
src/qseek/exporters/velest.py Show resolved Hide resolved
src/qseek/exporters/velest.py Outdated Show resolved Hide resolved
src/qseek/exporters/velest.py Outdated Show resolved Hide resolved
Copy link
Member

@miili miili left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Awesome! This is boiling down greatly :)

Comment on lines 67 to 68
with file.open("w") as fp:
fp.write(CONTROL_FILE_TPL.format(**self._asdict()))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use file.write_text(CONTROL_FILE_TPL.format(**self._asdict()) here.

Comment on lines 24 to 85
min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.3)

min_pick_semblance = FloatPrompt.ask("Minimum pick confidence", default=0.2)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rename to min_event_semblance. Use default values from self. Here self.min_event_semblance.

Comment on lines 177 to 180
if float(d) / 1000 not in dep_velest:
dep_velest.append(float(d) / 1000)
vp_velest.append(float(vp[i]) / 1000)
vs_velest.append(float(vs[i]) / 1000)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Instead of using 1000 introduce a constant KM = 1e3. This gives context what is happening here.

Consider using a zip: for d, vp, vs in zip(dep, vp, vs)

)
file.write("\n")
if count_p == 0 and count_s == 0:
logging.warning("Warning:No phases obesered for event{event.time}, removed")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Log suggestion:

logger.warning("Removing event {event.time}: No phases observed")

fpout.write("\n")

@staticmethod
def make_velmod_file(modname: Path, vp: list, vs: list, dep: list):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please complete typing for list... I suppose this is list[float]?

This can be a simple function, no need for a staticmethod.

Comment on lines 188 to 197
def velest_location(self, location: Location) -> str:
if location.effective_lat < 0:
velest_lat = f"{location.effective_lat:7.4f}S"
else:
velest_lat = f"{location.effective_lat:7.4f}N"
if location.effective_lon < 0:
velest_lon = f"{location.effective_lon:8.4f}W"
else:
velest_lon = f"{location.effective_lon:8.4f}E"
return velest_lat, velest_lon
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be a simple function, not a class method. self is not used.

Comment on lines 123 to 127
for event in catalog:
if event.semblance < min_pick_semblance:
continue
if event.receivers.n_observations(phase_p) < min_receivers_number:
continue
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also filter out events which are too close to the border.

Copy link
Member

@miili miili left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great, I have some last comments before merging.

Comment on lines 175 to 186
dep = search.ray_tracers.root[0].earthmodel.layered_model.profile("z")
vp = search.ray_tracers.root[0].earthmodel.layered_model.profile("vp")
vs = search.ray_tracers.root[0].earthmodel.layered_model.profile("vs")
dep_velest = []
vp_velest = []
vs_velest = []
km = 1000.0
for d, vpi, vsi in zip(dep, vp, vs, strict=False):
if float(d) / km not in dep_velest:
dep_velest.append(float(d) / km)
vp_velest.append(float(vpi) / km)
vs_velest.append(float(vsi) / km)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dont understand this part.

Why not:

depths = [depth/KM for depth in search.ray_tracers.root[0].earth
velocity_p = [vp/KM for vp in search.ray_tracers.root[0].earthmodel.layered_model.profile("vp")]
velocity_s = [vs/KM for vsin search.ray_tracers.root[0].earthmodel.layered_model.profile("vs")]

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi, thank you so much for the review! This is because the velocity model for velest requires only one value for each depth, but the velocity model of qseek is layered model with two velocity at one depth, so I need to remove another value.

Comment on lines 202 to 234
with outfile.open("a") as file:
file.write(
f"{event.time:%y%m%d %H%M %S.%f}"[:-4]
+ f" {lat} {lon} {event.depth/1000:7.2f} {mag:5.2f}\n"
)
count_p = 0
count_s = 0
for rec, dectection in observed_arrivals:
if dectection.observed.detection_value < 0.4:
quality_weight = 3
elif dectection.observed.detection_value < 0.6:
quality_weight = 2
elif dectection.observed.detection_value < 0.8:
quality_weight = 1
else:
quality_weight = 0
if dectection.phase.endswith("P"):
phase = "P"
count_p += 1
else:
phase = "S"
count_s += 1
traveltime = (dectection.observed.time - event.time).total_seconds()
file.write(
f" {rec.station:6s} {phase:1s} {quality_weight:1d} {traveltime:7.2f}\n"
)
file.write("\n")
if count_p == 0 and count_s == 0:
logger.warning("Removing event {event.time}: No phases observed")
with outfile.open("r") as file:
lines = file.readlines()
with outfile.open("w") as file:
file.writelines(lines[:-2])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here store the file output in a string

write_out = f"{event.time}..."
for rec_detection in observed_arrivals:
    write_out += f"..."

if count_p or count_s:
    with outfile.open("a") as file:
        file.write(write_out)     

allow_low_velocity: bool = False

def write_config_file(self, file: Path):
with file.open("w") as fp:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

open is not needed.

Comment on lines 87 to 101
min_receivers_number = IntPrompt.ask(
"Minimum number of receivers (P phase)", default=10
)
min_distance_to_border = FloatPrompt.ask(
"Minimum distance to border(m)", default=500.0
)
min_p_phase_confidence = FloatPrompt.ask(
"Minimum pick probability for P phase", default=0.3
)
min_s_phase_confidence = FloatPrompt.ask(
"Minimum pick probability for S phase", default=0.3
)
max_traveltime_delay = FloatPrompt.ask(
"Maximum difference between theoretical and observed arrival", default=2.5
)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use defaults from self

Comment on lines 210 to 217
if dectection.observed.detection_value < 0.4:
quality_weight = 3
elif dectection.observed.detection_value < 0.6:
quality_weight = 2
elif dectection.observed.detection_value < 0.8:
quality_weight = 1
else:
quality_weight = 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider using

np.digitize(dectection.observed.detection_value, [1.0, 0.8, 0.6, 0.4, 0.0])

@miili
Copy link
Member

miili commented Jul 16, 2024

Merged in dev

@miili miili closed this Jul 16, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants