diff --git a/.bumpversion.cfg b/.bumpversion.cfg index fd104aa..307bd1a 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,5 +1,5 @@ [bumpversion] -current_version = 1.2.1 +current_version = 1.3.0 commit = True tag = False parse = (?P\d+)\.(?P\d+)\.(?P\d+)(-(?P.*)-(?P\d+))? diff --git a/CHANGELOG.md b/CHANGELOG.md index a817b7b..c20dc00 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,10 @@ # Changelog + +## [1.3.0] 2025-12-04 +### Changed + - updated documentation on how to generate and plot MFD + - changed name of solvis module to utils and hoisted to top level namespace + ## [1.2.1] 2025-11-03 ### Changed - added `audit` to tox config diff --git a/docs/api/solvis/solvis.md b/docs/api/solvis/solvis.md deleted file mode 100644 index b7f38f3..0000000 --- a/docs/api/solvis/solvis.md +++ /dev/null @@ -1 +0,0 @@ -::: solvis.solvis \ No newline at end of file diff --git a/docs/api/solvis/utils.md b/docs/api/solvis/utils.md new file mode 100644 index 0000000..250b0f6 --- /dev/null +++ b/docs/api/solvis/utils.md @@ -0,0 +1 @@ +::: solvis.utils \ No newline at end of file diff --git a/docs/usage.md b/docs/usage.md index f4313d1..7e4e44f 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -106,40 +106,41 @@ filtered_rates_df.to_csv("filtered_rupture_rates.csv" ) Magnitude-Frequency Distribution (MFD) histogram data can also be generated and saved: -## Calculate an MFD with Pandas - +## Calculate and plot an MFD +This example will generate the incremental MFD from an inversion solution, calculate the cumulative MFD and plot them. ```py +import matplotlib.pyplot as plt +import numpy as np import pandas as pd -bins = [round(x / 100, 2) for x in range(500, 1000, 10)] -mfd = filtered_rates_df.groupby( - pd.cut( - filtered_rates_df.Magnitude, - bins=bins, - ))["rate_weighted_mean"].sum(numeric_only=False) -mfd.to_csv( "NSHM_filtered_mfd.csv") -``` -## Plot it using Matplotlib +from solvis.solution.inversion_solution import InversionSolution +from solvis.utils import mfd_hist -```py -import matplotlib.pyplot as plt -import numpy as np +solution = InversionSolution.from_archive("InversionSolution.zip") +mfd = mfd_hist(solution.model.ruptures_with_rupture_rates) + + +def plot_mfd(mfd: pd.DataFrame, ax): -def plot_mfd(mfd: pd.DataFrame, title: str = "Title"): + full_index = pd.CategoricalIndex(mfd.index.categories) + full_mfd = pd.Series(index=full_index, data=0) + mfd_cumulative = full_mfd + mfd + mfd_cumulative.fillna(0, inplace=True) + + mfd_cumulative = mfd_cumulative.loc[::-1].cumsum().loc[::-1] mag = [a.mid for a in mfd.index] + mag_cumulative = [a.mid for a in mfd_cumulative.index] + dmag = mag[1] - mag[0] rate = np.asarray(mfd) rate[rate == 0] = 1e-20 # set minimum rate for log plots - fig = plt.figure() - fig.set_facecolor("white") - plt.title(title) - plt.ylabel("Incremental Rate ($yr^-1$)") - plt.xlabel("Magnitude") - plt.semilogy(mag, rate, color="red") - plt.axis([6.0, 9.0, 0.000001, 1.0]) - plt.grid(True) - return plt - -plot = plot_mfd(mfd, title="A Filtered MFD") -plot.savefig("A filtered MFD plot.png") -plot.close() + ax.set_yscale('log') + ax.bar(x=mag, height=rate, width=0.8 * dmag) + ax.plot(mag_cumulative, mfd_cumulative) + ax.set_ylim([1e-10, 1]) + + +fig = plt.figure() +ax = plt.gca() +plot_mfd(mfd, ax) +plt.show() ``` \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index e080b3c..d84a259 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -38,7 +38,7 @@ nav: - typing: api/solution/typing.md - config: api/solvis/config.md - geometry: api/solvis/geometry.md - - solvis: api/solvis/solvis.md + - utils: api/solvis/utils.md - Development and Contributing: - contributing.md - Testing: testing.md diff --git a/pyproject.toml b/pyproject.toml index a4ab6a2..da07454 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "solvis" -version = "1.2.1" +version = "1.3.0" description = "analysis of opensha modular solution files." readme = "README.md" license = "AGPL-3.0-or-later" diff --git a/solvis/__init__.py b/solvis/__init__.py index 9628b54..8a97b41 100644 --- a/solvis/__init__.py +++ b/solvis/__init__.py @@ -10,6 +10,7 @@ CompositeSolution: the container class for complete model and logic tree. """ +from . import utils from .solution import CompositeSolution, FaultSystemSolution, InversionSolution -__version__ = '1.2.1' +__version__ = '1.3.0' diff --git a/solvis/scripts/cli.py b/solvis/scripts/cli.py index ea0b577..e2e6d69 100644 --- a/solvis/scripts/cli.py +++ b/solvis/scripts/cli.py @@ -16,7 +16,7 @@ from solvis.geometry import circle_polygon from solvis.get_secret import get_secret from solvis.solution.inversion_solution import BranchInversionSolution, InversionSolution -from solvis.solvis import export_geojson +from solvis.utils import export_geojson # Get API key from AWS secrets manager API_URL = os.getenv('NZSHM22_TOSHI_API_URL', "http://127.0.0.1:5000/graphql") diff --git a/solvis/solvis.py b/solvis/solvis.py deleted file mode 100644 index 923d4d8..0000000 --- a/solvis/solvis.py +++ /dev/null @@ -1,26 +0,0 @@ -#!python3 -""" -Some original solvis API helper functions are defined in this module. - -But most are now deprecated and replaced by equivalent filter and model methods. -""" -from pathlib import Path -from typing import Union - -import geopandas as gpd -import pandas as pd - - -def mfd_hist(ruptures_df: pd.DataFrame, rate_column: str = "Annual Rate"): - # https://stackoverflow.com/questions/45273731/binning-a-column-with-python-pandas - bins = [round(x / 100, 2) for x in range(500, 1000, 10)] - # Added observed=True in advance of default change (from False) as advised in pandas FutureWarning - mfd = ruptures_df.groupby(pd.cut(ruptures_df.Magnitude, bins=bins), observed=True)[rate_column].sum() - return mfd - - -def export_geojson(gdf: gpd.GeoDataFrame, filename: Union[str, Path], **kwargs): - print(f"Exporting to {filename}") - f = open(filename, 'w') - f.write(gdf.to_json(**kwargs)) - f.close() diff --git a/solvis/utils.py b/solvis/utils.py new file mode 100644 index 0000000..27a27a5 --- /dev/null +++ b/solvis/utils.py @@ -0,0 +1,43 @@ +"""Helper functions for solvis.""" +from pathlib import Path +from typing import Union + +import geopandas as gpd +import pandas as pd + + +def mfd_hist(ruptures_df: pd.DataFrame, rate_column: str = "Annual Rate") -> pd.Series: + """Generate a magnitude frequency distribution (MFD) from a DataFrame of ruptures. + + Args: + ruptures_df: dataframe of ruptures + rate_column: the dataframe column name containing the ruputre rates + + Returns: + a pandas Series of the MFD histogram. + + Example: + ```py + solution = InversionSolution.from_archive("InversionSolution.zip") + mfd = mfd_hist(solution.model.ruptures_with_rupture_rates) + ``` + """ + # https://stackoverflow.com/questions/45273731/binning-a-column-with-python-pandas + bins = [round(x / 100, 2) for x in range(500, 1000, 10)] + # Added observed=True in advance of default change (from False) as advised in pandas FutureWarning + mfd = ruptures_df.groupby(pd.cut(ruptures_df.Magnitude, bins=bins), observed=True)[rate_column].sum() + return mfd + + +def export_geojson(gdf: gpd.GeoDataFrame, filename: Union[str, Path], **kwargs): + """Export GeoDataFrame to json file. + + Args: + gdf: a GeoDataFrame to export to file + filename: the path of the file: + **kwargs: key-word arguments to pass to the GeoDataFrame.to_json() function. + """ + print(f"Exporting to {filename}") + f = open(filename, 'w') + f.write(gdf.to_json(**kwargs)) + f.close() diff --git a/test/test_solvis_functions.py b/test/test_solvis_functions.py index 9e3461a..b72d048 100644 --- a/test/test_solvis_functions.py +++ b/test/test_solvis_functions.py @@ -1,13 +1,12 @@ import pandas as pd import pytest -import solvis.solvis +import solvis.utils def test_mfd_hist(crustal_small_fss_fixture, crustal_solution_fixture): - - mfd = solvis.solvis.mfd_hist(crustal_small_fss_fixture.model.ruptures_with_rupture_rates, "rate_weighted_mean") + mfd = solvis.utils.mfd_hist(crustal_small_fss_fixture.model.ruptures_with_rupture_rates, "rate_weighted_mean") assert mfd.loc[pd.Interval(left=7.0, right=7.1)] == pytest.approx(0.0011956305) - mfd = solvis.solvis.mfd_hist(crustal_solution_fixture.model.ruptures_with_rupture_rates) + mfd = solvis.utils.mfd_hist(crustal_solution_fixture.model.ruptures_with_rupture_rates) assert mfd.loc[pd.Interval(left=7.1, right=7.2)] == pytest.approx(0.0018980678)