Skip to content

Alpine plots#107

Draft
lispandfound wants to merge 41 commits intomasterfrom
alpine_plots
Draft

Alpine plots#107
lispandfound wants to merge 41 commits intomasterfrom
alpine_plots

Conversation

@lispandfound
Copy link
Contributor

No description provided.

@gemini-code-assist
Copy link
Contributor

Summary of Changes

Hello @lispandfound, I'm Gemini Code Assist1! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request introduces a comprehensive suite of new plotting tools and significantly enhances existing visualization capabilities. The primary goal is to provide more robust and versatile methods for analyzing and presenting seismic simulation outputs, particularly in comparison with empirical Ground Motion Models. Key improvements include advanced 3D ground motion animations, detailed GMM comparison plots, and the integration of R for statistical data fitting. These changes streamline the process of generating high-quality visualizations and facilitate deeper scientific analysis of simulation results.

Highlights

  • New Plotting Tools: Introduced several new scripts for visualizing seismic data, including Ground Motion Model (GMM) comparisons, pSA vs. Rrup plots, Fourier Amplitude Spectra (FAS) plots, response spectra, and raw waveforms. These new tools significantly expand the project's analytical and presentation capabilities.
  • 3D Ground Motion Animation with PyVista: The plot_ts script has been completely refactored to utilize pyvista for advanced 3D visualization of ground motion. This replaces the previous 2D matplotlib/ffmpeg approach, offering a more immersive and detailed animated experience, including Digital Elevation Model (DEM) integration and dynamic fault plane rendering.
  • GMM Integration and Analysis: New utilities and scripts have been added to compute and plot comparisons between simulation results and NSHM2022 Ground Motion Model predictions. This includes functions to handle site and rupture properties, enabling deeper insights into the agreement between simulations and empirical models.
  • R Integration for LOESS Fitting: The rpy2 library has been added as a new dependency, enabling seamless integration with R for advanced statistical methods. Specifically, LOESS (Locally Estimated Scatterplot Smoothing) fitting is now available for data analysis within plotting scripts, enhancing the ability to identify trends in noisy data.
  • Enhanced Plotting Customization and Utilities: Existing plotting functions (plot_realisation, plot_srf, plot_sources) have been refined for better control over map features, source representation (including fault traces), and general plot aesthetics. A new balanced_subplot_grid utility simplifies the creation of flexible subplot layouts.
Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature Command Description
Code Review /gemini review Performs a code review for the current pull request in its current state.
Pull Request Summary /gemini summary Provides a summary of the current pull request in its current state.
Comment @gemini-code-assist Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help /gemini help Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Footnotes

  1. Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution.

Copy link
Contributor

@gemini-code-assist gemini-code-assist bot left a comment

Choose a reason for hiding this comment

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

Code Review

This pull request introduces several new plotting scripts for visualizing simulation results, including GMM comparisons, response spectra, FAS, and waveforms. It also refactors existing plotting scripts to improve code quality and add new features. Review comments highlight several issues, including hardcoded paths and values, type mismatches, incorrect calculations, incomplete documentation, and deprecated Matplotlib arguments. The reviewer suggests making paths and parameters configurable, correcting type mismatches, using named constants for magic numbers, completing docstrings, and ensuring compatibility with older Python versions.

lines.extend([point_a, point_b])

cmap_min, cmap_max, cmap = cmap_from_cpt(
Path("/home/jake/tmp/palm_springs_nz_topo.cpt")
Copy link
Contributor

Choose a reason for hiding this comment

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

critical

The path "/home/jake/tmp/palm_springs_nz_topo.cpt" is hardcoded and points to a specific user's temporary directory. This makes the code non-portable and will cause failures on other systems. This path must be made configurable, either as an argument or a setting, and should ideally point to a resource within the project or a well-defined location.

    # TODO: Make CPT file path configurable
    cmap_min, cmap_max, cmap = cmap_from_cpt(
        Path("/path/to/your/cpt/file.cpt")
    )

) -> None:
freqs = dataset.frequency.values
vs30 = dataset.vs30.sel(station=station).item()
site_properties = utils.compute_site_properties(vs30)
Copy link
Contributor

Choose a reason for hiding this comment

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

high

The utils.compute_site_properties function expects an xr.Dataset but is being called with a single vs30 float. This is a type mismatch and will cause a runtime error. Please adjust the call to pass the entire site xarray dataset, or modify utils.compute_site_properties to accept a float if that's the intended usage here.

Suggested change
site_properties = utils.compute_site_properties(vs30)
site_properties = utils.compute_site_properties(site)

avg_rake = np.degrees(circmean(np.radians(rakes), np.array(moments)))
avg_dip = float(np.average(dips, weights=moments))
avg_moment = float(np.mean(moments))
total_moment = avg_moment * len(moments)
Copy link
Contributor

Choose a reason for hiding this comment

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

high

The calculation total_moment = avg_moment * len(moments) is incorrect. If avg_moment is already the mean of the individual moments, then total_moment should be the sum of the moments, which is np.sum(moments). Multiplying the average by the count of moments will give the sum, but it's more direct and clearer to use np.sum if moments is a list/array of individual moments.

Suggested change
total_moment = avg_moment * len(moments)
total_moment = np.sum(moments)

Comment on lines +418 to +436
def compute_site_properties(sites: xr.Dataset) -> SiteProperties:
vs30 = sites.vs30.values
z1pt0 = oqw.estimations.chiou_young_08_calc_z1p0(vs30)
z2pt5 = oqw.estimations.chiou_young_08_calc_z2p5(vs30)
rrup = sites.rrup.values
rjb = sites.rjb.values
rx = rjb
ry = rjb
vs30measured = False
return SiteProperties(
vs30measured=vs30measured,
vs30=vs30,
z1pt0=z1pt0,
z2pt5=z2pt5,
rrup=rrup,
rjb=rjb,
rx=rx,
ry=ry,
)
Copy link
Contributor

Choose a reason for hiding this comment

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

high

The compute_site_properties function is designed to take an xr.Dataset as input, but it appears to be called with a single vs30 float in visualisation/waveforms/plot_fas.py. This is a type mismatch that will lead to a runtime error. The function either needs to be refactored to accept a single vs30 value, or the calling code needs to pass an xr.Dataset containing the site properties.

def compute_site_properties(
    vs30: float | npt.NDArray[np.floating],
    rrup: float | npt.NDArray[np.floating],
    rjb: float | npt.NDArray[np.floating],
) -> SiteProperties:
    vs30_val = vs30 if isinstance(vs30, (float, np.floating)) else vs30.item()
    z1pt0 = oqw.estimations.chiou_young_08_calc_z1p0(vs30_val)
    z2pt5 = oqw.estimations.chiou_young_08_calc_z2p5(vs30_val)
    rx = rjb
    ry = rjb
    vs30measured = False
    return SiteProperties(
        vs30measured=vs30measured,
        vs30=vs30,
        z1pt0=z1pt0,
        z2pt5=z2pt5,
        rrup=rrup,
        rjb=rjb,
        rx=rx,
        ry=ry,
    )

# Write the raw video data to FFmpeg's stdin
for frame in frames:
process.stdin.write(frame)
scalars[scalars < 10] = np.nan
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The value 10 is a magic number. Please define it as a named constant with a clear explanation of its purpose (e.g., MIN_SCALAR_THRESHOLD) to improve readability and maintainability.

        MIN_SCALAR_THRESHOLD = 10
        scalars[scalars < MIN_SCALAR_THRESHOLD] = np.nan

ax.set_xlim(freqs.min(), 50)
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(visible=True, which="both", axis="both", lw=0.3)
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The visible argument for ax.grid() is deprecated in Matplotlib 3.4 and later. Please remove visible=True as True is the default behavior for ax.grid().

Suggested change
ax.grid(visible=True, which="both", axis="both", lw=0.3)
ax.grid(which="both", axis="both", lw=0.3)



def find_region(domain: DomainParameters) -> tuple[float, float, float, float]:
"""Find an appropriate domain,"""
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The docstring for find_region is incomplete. It ends with a comma, suggesting that the description is unfinished. Please complete the docstring to clearly explain what the function does.

Suggested change
"""Find an appropriate domain,"""
def find_region(domain: DomainParameters) -> tuple[float, float, float, float]:
"""Find an appropriate domain for plotting."""

/ 1000
)
rupture_df = pd.DataFrame(
{
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The value vs30measured=False is hardcoded. If this value can vary, it should be derived from the input data or made a configurable parameter. If it's always False for this context, a comment explaining why would be beneficial.

            "vs30measured": False, # Assuming vs30 is not measured for this context

ax.set_xlabel("Period [s]")
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(visible=True, which="both", axis="both", lw=0.3)
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The visible argument for ax.grid() is deprecated in Matplotlib 3.4 and later. Please remove visible=True as True is the default behavior for ax.grid().

Suggested change
ax.grid(visible=True, which="both", axis="both", lw=0.3)
ax.grid(which="both", axis="both", lw=0.3)

CMS2 = "cm/s^2"


_UNIT_CONVERSION_TABLE = {Units.G: 981.0, Units.CMS: 1.0, Units.CMS2: 1.0}
Copy link
Contributor

Choose a reason for hiding this comment

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

medium

The conversion factor 981.0 for Units.G is hardcoded. While this is the standard acceleration due to gravity in cm/s², it's a magic number. Consider defining it as a named constant (e.g., GRAVITY_CMS2) to improve readability and maintainability.

Suggested change
_UNIT_CONVERSION_TABLE = {Units.G: 981.0, Units.CMS: 1.0, Units.CMS2: 1.0}
GRAVITY_CMS2 = 981.0
_UNIT_CONVERSION_TABLE = {Units.G: GRAVITY_CMS2, Units.CMS: 1.0, Units.CMS2: 1.0}

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.

1 participant