107 mfd hist

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.2.1
+current_version = 1.3.0
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+)\.(?P<patch>\d+)(-(?P<release>.*)-(?P<build>\d+))?

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

docs/api/solvis/utils.md

@@ -0,0 +1 @@
+::: solvis.utils

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()
 ```

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

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"

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'

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")

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()

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)