Merging Wenjie's edits from years ago #20
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
…od and data_fit method
Two tests are still not passed: 1) code formatting: going to rebase to fix this problem 2) figure for tests should be replaced in the future
Conflicts: src/pyflex/tests/test_pyflex.py src/pyflex/window_selector.py
…plotlib(should be added back later)
…tion 1) modify `calculate_preliminary` to check the signal region instead 2) add `reject_on_noise_region`, split from `reject_on_selection_mode` 3) add selection_mode=`custom` for user customization
Trying to merge the merge branch with Wenjie's edit into the main branch
|
Thanks @lsawade! Good to have those changes merged back into the repo, I appreciate you making the effort! |
|
Hey @lsawade, just wanted to let you know this is still on my to-do list! Apologies for the long wait |
|
All good! I found a bug that was introduced by the merged anyways! |
|
Accidentally clicked resolve... |
|
Have to do it tonight, traveling rn |
|
Woops, sorry if you're getting emails for all these comments, I'm just going through the PR now, will leave a longer message when I finish looking through. And no rush! edit: that was my bad, I was leaving single comments and not putting all into one review, hopefully you do not get inundated by future stuff! |
bch0w
requested changes
Jan 29, 2024
There was a problem hiding this comment.
Hi @lsawade, thanks again for pushing these through and sorry again for how long it's taken me to get to it! Also sorry for harassing you with notifications when I first started the review :D
I've just gone through the diffs and looks great, I've left code change suggestions on some things, but otherwise this looks like this is a pretty useful and non-intrusive PR. After those changes I am happy to merge through and version push to 0.3.0.
Just summarizing the changes I saw so I can copy paste this into the release notes:
- Adds new Config parameters:
s2n_limit_energy,max_surface_wave_velocity,max_time_after_last_arrival,selection_mode,window_signal_to_noise_type - Adds ability to reject windows based on signal to noise energy ratio (original code based on amplitudes)
- Feature: force window selection to certain phase arrivals with
selection_mode - Adds new Window parameters:
channel_id_2,snr_amplitude,snr_amplitude_threshold,distance_in_deg - WindowSelector.write_text provides a simple plaintext file write function
- Adds a new rejection function to reject peaks and troughs not within the desired signal region which is determined by
selection_modeand expected arrival times - Adds new debug log message to list out chosen windows
- Adds a new rejection criteria for windows in pre-defined noise regions
Let me know what you think, and no rush on this. Thanks again!
| :type min_surface_wave_velocity: float | ||
|
|
||
| :param max_surface_wave_velocity: The maxium surface wave velocity |
There was a problem hiding this comment.
Suggested change
| :param max_surface_wave_velocity: The maxium surface wave velocity | |
| :param max_surface_wave_velocity: The maximum surface wave velocity |
| "6) custom; 7) phase_list;".format(selection_mode)) | ||
|
|
||
| self.selection_mode = selection_mode | ||
|
|
Comment on lines
+84
to
+85
|
|
||
| Returns the list of surface wave arrivals in time order. |
There was a problem hiding this comment.
Suggested change
| Returns the list of surface wave arrivals in time order. | |
| :type dist_in_deg: float | |
| :param dist_in_deg: great circle source-receiver distance in units of degrees | |
| :type min_vel: float | |
| :param min_vel: minimum expected surface wave velocity to get latest arrival, units of km/s | |
| :type max_vel: float | |
| :param max_vel: maximum expected surface wave velocity to get earliest arrival, units of km/s | |
| :type ncircles: int | |
| :param ncircles: number of global orbits to calculate, defaults to 1 | |
| Returns the list of surface wave arrivals in time order. |
| Calculate the arrival time of surface waves, based on the distance | ||
| and velocity range (min_vel, max_vel). | ||
| This function will calculate both minor-arc and major-arc surface | ||
| waves. It further calcualte the surface orbit multiple times |
There was a problem hiding this comment.
Suggested change
| waves. It further calcualte the surface orbit multiple times | |
| waves. It further calculates the surface orbit multiple times |
| Returns the list of surface wave arrivals in time order. | ||
| """ | ||
| if min_vel > max_vel: | ||
| min_vel, max_vel = max_vel, min_vel |
There was a problem hiding this comment.
Maybe we should add a logger warning message to let the User know they mixed these up?
Comment on lines
+765
to
+849
| def _prepare_min_max_timebox(self, min_time, max_time, srate): | ||
| logger.debug("selection_timebox min/max time: {:.2f} {:.2f}" | ||
| "".format(min_time, max_time)) | ||
|
|
||
| timebox = np.zeros(self.observed.stats.npts, dtype=bool) | ||
| il = int(min_time * srate) | ||
| ir = int(max_time * srate) | ||
| timebox[il:(ir+1)] = 1 | ||
|
|
||
| logger.debug("selection_timebox index range: {} - {} / {}".format( | ||
| il, ir, len(timebox))) | ||
|
|
||
| self.selection_timebox = timebox | ||
| return timebox | ||
|
|
||
| def _prepare_exclude_surface_wave_timebox(self, min_time, max_time, | ||
| srate, offset): | ||
| logger.debug("selection_timebox exclude surface waves") | ||
|
|
||
| # get the surface windows first | ||
| self._prepare_surface_wave_timebox(srate, offset) | ||
|
|
||
| # flip the sign to reject surface windoes | ||
| timebox = self.selection_timebox | ||
| timebox = ~timebox | ||
|
|
||
| # set the min_time and max_time range to non-selectable | ||
| i1 = int(min_time * srate) | ||
| timebox[:i1] = False | ||
|
|
||
| i2 = int(max_time * srate) | ||
| timebox[i2:] = False | ||
|
|
||
| logger.debug("Selectable region coverage percentage: {}/{}".format( | ||
| timebox.sum(), len(timebox))) | ||
| self.selection_timebox = timebox | ||
| return timebox | ||
|
|
||
| def _prepare_surface_wave_timebox(self, srate, offset): | ||
| logger.debug("selection_timebox using surface waves") | ||
|
|
||
| max_vel = self.config.max_surface_wave_velocity | ||
| min_vel = self.config.min_surface_wave_velocity | ||
| arrivals = get_surface_wave_arrivals( | ||
| self.dist_in_deg, min_vel, max_vel, ncircles=2) | ||
| logger.debug("Number of surface wave arrivals: {}".format( | ||
| len(arrivals))) | ||
|
|
||
| timebox = np.zeros(self.observed.stats.npts, dtype=bool) | ||
|
|
||
| for _i, arr in enumerate(arrivals): | ||
| il = int((arr[0] + offset - self.config.min_period) * srate) | ||
| ir = int((arr[1] + offset + 2.5 * self.config.min_period) * srate) | ||
| logger.debug(f"Arrival {_i} (min, max): ({il}, {ir})") | ||
| timebox[il:(ir+1)] = 1 | ||
| logger.debug("Selectable region coverage percentage: {}/{}".format( | ||
| timebox.sum(), len(timebox))) | ||
| self.selection_timebox = timebox | ||
| return timebox | ||
|
|
||
| def _prepare_phase_timebox(self, phase_list, srate, buffer_time, offset): | ||
| logger.debug("Selection_timebox using phases: {}".format(phase_list)) | ||
|
|
||
| source_depth = self.event.depth_in_m / 1000.0 | ||
| logger.debug("event depth: {:.2f} km".format(source_depth)) | ||
| logger.debug("distance in deg: {:.2f}".format(self.dist_in_deg)) | ||
| tts = self.taupy_model.get_travel_times( | ||
| source_depth_in_km=source_depth, | ||
| distance_in_degree=self.dist_in_deg, | ||
| phase_list=phase_list | ||
| ) | ||
|
|
||
| timebox = np.zeros(self.observed.stats.npts, dtype=bool) | ||
| for t in tts: | ||
| t_left = t.time - buffer_time + offset | ||
| il = int(t_left * srate) | ||
| t_right = t.time + buffer_time + offset | ||
| ir = int(t_right * srate) | ||
| timebox[il:(ir+1)] = 1 | ||
|
|
||
| logger.debug("Selectable region percentage: {}/{}".format( | ||
| timebox.sum(), len(timebox))) | ||
| self.selection_timebox = timebox | ||
|
|
||
| return timebox |
There was a problem hiding this comment.
I think this can be condensed, alongside reject_on_selection_criteria. I think it's fine for now but in the future I/we might want to try and find a cleaner way to write this.
Comment on lines
+877
to
+881
| Reject based on selection mode. | ||
| This function will reject windows outside of the | ||
| wave category specified. For example, if config.selection_mode | ||
| == "body_waves", only body wave windows will be selected(after | ||
| first arrival and before surface wave arrival). |
There was a problem hiding this comment.
Reject windows outside a pre-defined phase window (e.g., body waves,
surface waves, mantle waves) using the User-defined `selection_mode`.
Availalable values for `selection_mode` are:
- all_waves: all phases after first arrival
- body_waves: all phases after first arrival and before first surface
- surface_waves: all phases inside surface wave train
- mantle_waves: all phases after final surface wave arrival
- body_and_mantle_waves: exclude surface waves
- body_and_surface_waves: all phases before final surface wave arrival
Comment on lines
+260
to
+261
| * ``"body_and_surface_waves"``: windows inside body and surface | ||
| wave regions. |
There was a problem hiding this comment.
Suggested change
| * ``"body_and_surface_waves"``: windows inside body and surface | |
| wave regions. | |
| * ``"body_and_surface_waves"``: windows inside body and surface | |
| wave regions. | |
| * ``"body_and_mantle_waves"``: exclude surface waves |
| win.cc_shift) | ||
| if not (tshift_min <= win.cc_shift_in_seconds <= tshift_max): | ||
| logger.debug("Window [%.1f - %.1f] rejected due to time " | ||
| "shift does not satisfy: %.1f < %.1f < %.1f" |
There was a problem hiding this comment.
This is great, I've been meaning to change that log message because it's getting annoying having to manually convert samples to time every inversion!
| @@ -870,7 +1294,7 @@ def _sanity_checks(self): | |||
|
|
|||
| ptp = sorted([self.observed.data.ptp(), self.synthetic.data.ptp()]) | |||
| if ptp[1] / ptp[0] >= 5: | |||
| warnings.warn("The amplitude difference between data and " | |||
| warnings.warn("The amplitude difference between data and" | |||
There was a problem hiding this comment.
Suggested change
| warnings.warn("The amplitude difference between data and" | |
| warnings.warn("The amplitude difference between data and " |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Hi @bch0w !
I'm reorganizing some of my packages, one of which depends on some of the edits Wenjie made years ago. I have merged them in a hopefully non-conflicting way. They work with my workflows, and tests are passing, but I think it would be good if you double check this with your workflows.
Cheers,
Lucas