WIP: Fix sampling, other mbt

openquake/cat/completeness/mfd_eval_plots.py

@@ -23,6 +23,7 @@ def _read_results(results_dir):
         dfs.append(df)
 
     df_all = pd.concat(dfs, ignore_index=True)
+    df_all = df_all[df_all['agr'].notna()]
 
     mags, rates = [], []
     cm_rates = []
@@ -168,7 +169,10 @@ def plot_best_mfds(df_best, figsdir):
             if num <= 10:
                 alpha1 = 0.1
             else:
-                alpha1 = 10/num
+                alpha1 = 2/num
+
+            if alpha1 > 0.2:
+                breakpoint()
 
             plt.scatter(row.mags, row.rates, marker='_', color='r', 
                         alpha=alpha1)
@@ -219,7 +223,7 @@ def plot_weighted_covariance_ellipse(df, figdir, n_std=1.0,
     weights = (wei - min(wei)) / (max(wei) - min(wei)) 
 
     # set up plot
-    fig, ax = plt.subplots()
+    fig, ax = plt.subplots(figsize=(10,10))
     hb = ax.hexbin(x, y, gridsize=20, cmap='Blues')
     cb = fig.colorbar(hb, ax=ax)
     cb.set_label('counts')
@@ -276,9 +280,11 @@ def plot_weighted_covariance_ellipse(df, figdir, n_std=1.0,
             fontsize=12, color=color)
 
 
-    ax.set_xlabel('a-value', fontsize=12)
-    ax.set_ylabel('b-value', fontsize=12)
-    ax.set_title(figname.replace('.png', ''))
+    ax.set_xlabel('a-value', fontsize=16)
+    ax.set_ylabel('b-value', fontsize=16)
+    ax.set_title(figname.replace('.png', ''), fontsize=16)
+    ax.xaxis.set_tick_params(labelsize=14)
+    ax.yaxis.set_tick_params(labelsize=14)
 
     fout = os.path.join(figdir, figname)
     plt.savefig(fout, dpi=300)
@@ -290,39 +296,41 @@ def plot_weighted_covariance_ellipse(df, figdir, n_std=1.0,
 
 def plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=10, bw=0.2, figname=None):
 # Group the DataFrame by the 'Category' column and apply the histogram calculation function
+
+    fig, ax = plt.subplots(figsize=(10,6))
 
     fl_mags = [item for sublist in df_all.mags.values for item in sublist]
     fl_rates = [item for sublist in df_all.rates.values for item in sublist]
     fl_crates = [item for sublist in df_all.cm_rates.values for item in sublist]
     fl_df = pd.DataFrame({'mags': fl_mags, 'rates': fl_rates, 'cm_rates': fl_crates})
-
+     
     grouped = fl_df.groupby('mags')
     hist_data = grouped.apply(lambda g: norm_histo(g, field=field, bins=bins))
     mags = hist_data._mgr.axes[0].values
-    results = hist_data.values
+    #results = hist_data.values
 
-    all_alpha = []
-    for mag, rat in zip(mags, results):
-        m = [mag] * len(rat[0])
-        nmc = rat[1]
-        all_alpha.extend(nmc)
-    alph_min = min(all_alpha)
-    alph_max = max(all_alpha)
+    #all_alpha = []
+    #for mag, rat in zip(mags, results):
+    #    m = [mag] * len(rat[0])
+    #    nmc = rat[1]
+    #    all_alpha.extend(nmc)
+    #alph_min = min(all_alpha)
+    #alph_max = max(all_alpha)
 
-    norm = mcolors.Normalize(vmin=alph_min, vmax=alph_max)
+    #norm = mcolors.Normalize(vmin=alph_min, vmax=alph_max)
 
     # Choose a colormap
-    colormap = plt.cm.Purples
+    #colormap = plt.cm.Purples
 
-    for mag, rat in zip(mags, results):
-        m = [mag] * len(rat[0])
-        alph = rat[1] 
-        alph[alph<0.1]=0.2
-        plt.semilogy([m[0], m[0]], [min(rat[0]), max(rat[0])], c='gray', linewidth=1, zorder=1)
-        plt.scatter(m, rat[0], 250, marker='_', c=alph, cmap=colormap, norm=norm, zorder=0)
+    #for mag, rat in zip(mags, results):
+    #    m = [mag] * len(rat[0])
+    #    alph = rat[1] 
+    #    alph[alph<0.1]=0.2
+    #    ax.semilogy([m[0], m[0]], [min(rat[0]), max(rat[0])], c='gray', linewidth=1, zorder=1)
+    #    ax.scatter(m, rat[0], 250, marker='_', c=alph, cmap=colormap, norm=norm, zorder=0)
 
     for index, row in df_best.iterrows():
-        plt.scatter(row['mags'], row[field], 2, 'k', marker='s')
+        ax.scatter(row['mags'], row[field], 2, 'k', marker='s')
         mfd = TruncatedGRMFD(min(mags)-bw, 8.5, bw, row.agr, row.bgr)
         mgrts = mfd.get_annual_occurrence_rates()
         mfd_m = [m[0] for m in mgrts]
@@ -334,20 +342,22 @@ def plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=10, bw=0.2, fign
 
         if field == 'cm_rates':
             mfd_cr = [sum(mfd_r[ii:]) for ii in range(len(mfd_r))]
-            plt.semilogy(mfd_m, mfd_cr, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)
+            ax.semilogy(mfd_m, mfd_cr, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)
         else:
-            plt.semilogy(mfd_m, mfd_r, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)
+            ax.semilogy(mfd_m, mfd_r, color='maroon', linewidth=0.2, zorder=10, alpha=alpha)
 
     if figname==None:
         figname = f'mfds_all_{field}.png'
 
     fout = os.path.join(figsdir, figname)
 
 
-    plt.xlabel('Magnitude')
-    plt.ylabel('annual occurrence rate')
-    plt.title(figname.replace('.png', ''))
-    plt.savefig(fout)
+    ax.set_xlabel('Magnitude', fontsize=16)
+    ax.set_ylabel('annual occurrence rate', fontsize=16)
+    ax.set_title(figname.replace('.png', ''), fontsize=16)
+    ax.xaxis.set_tick_params(labelsize=14)
+    ax.yaxis.set_tick_params(labelsize=14)
+    plt.savefig(fout, dpi=300)
     plt.close()
 
 
@@ -376,16 +386,17 @@ def make_all_plots(resdir_base, compdir, figsdir_base, labels):
                         weights = nm_weight)
         print('plotting mfds')
         plot_best_mfds(df_best, figsdir)
+
         plot_all_mfds(df_all, df_best, figsdir, field='rates', bins=60)
         plot_all_mfds(df_all, df_best, figsdir, field='cm_rates', bins=60)
         plot_all_mfds(df_best, df_best, figsdir, field='rates', bins=60, 
                       figname='mfds_best_rates.png')
         plot_all_mfds(df_best, df_best, figsdir, field='cm_rates', bins=60, 
                       figname='mfds_best_cmrates.png')
-        plot_all_mfds(df_all, df_thresh, figsdir, field='rates', bins=60, 
-                      figname='mfds_thresh_rates.png')
-        plot_all_mfds(df_all, df_thresh, figsdir, field='cm_rates', bins=60, 
-                      figname='mfds_thresh_cmrates.png')
+        #plot_all_mfds(df_all, df_thresh, figsdir, field='rates', bins=60, 
+        #              figname='mfds_thresh_rates.png')
+        # plot_all_mfds(df_all, df_thresh, figsdir, field='cm_rates', bins=60, 
+        #              figname='mfds_thresh_cmrates.png')
         print('plotting covariance')
         cx, cy, mx1, my1, mx2, my2 = plot_weighted_covariance_ellipse(df_best, figsdir)
         plot_weighted_covariance_ellipse(df_thresh, figsdir, figname=f'{label}-20percent.png')

openquake/fnm/fault_modeler.py

@@ -340,10 +340,13 @@ def subdivide_rupture_mesh(
             subsec_lats = lats[i_start:i_end, j_start:j_end]
             subsec_depths = depths[i_start:i_end, j_start:j_end]
 
-            subsec_mesh = RectangularMesh(
-                subsec_lons, subsec_lats, subsec_depths
-            )
-            subsec_meshes.append({'row': i, 'col': j, 'mesh': subsec_mesh})
+            try:
+                subsec_mesh = RectangularMesh(
+                    subsec_lons, subsec_lats, subsec_depths
+                )
+                subsec_meshes.append({'row': i, 'col': j, 'mesh': subsec_mesh})
+            except:
+                print(i_start, i_end, j_start, j_end, i, j)
 
             j_start += n_subsec_pts_strike - 1
         i_start += n_subsec_pts_dip - 1
@@ -697,6 +700,7 @@ def get_boundary_3d(smsh):
         )
         for i in idx
     ]
+    tmp = [c for c in tmp if c != (0.0, 0.0, -0.0)]
     trace = LineString(tmp)
     coo.extend(tmp)
 
@@ -710,6 +714,7 @@ def get_boundary_3d(smsh):
         )
         for i in idx
     ]
+    tmp = [c for c in tmp if c != (0.0, 0.0, -0.0)]
     coo.extend(tmp)
 
     # Lower boundary
@@ -722,6 +727,7 @@ def get_boundary_3d(smsh):
         )
         for i in np.flip(idx)
     ]
+    tmp = [c for c in tmp if c != (0.0, 0.0, -0.0)]
     coo.extend(tmp)
 
     # Left boundary
@@ -734,6 +740,7 @@ def get_boundary_3d(smsh):
         )
         for i in np.flip(idx)
     ]
+    tmp = [c for c in tmp if c != (0.0, 0.0, -0.0)]
     coo.extend(tmp)
 
     return trace, Polygon(coo)
@@ -758,7 +765,10 @@ def make_subfault_gdf(subfault_df, keep_surface=False, keep_trace=False):
 
 
 def make_rupture_gdf(
-    fault_network, rup_df_key='rupture_df', keep_sequences=False
+    fault_network,
+    rup_df_key='rupture_df',
+    keep_sequences=False,
+    same_size_arrays: bool = True,
 ) -> gpd.GeoDataFrame:
     """
     Makes a GeoDataFrame, with a row for each rupture in the fault network.
@@ -785,7 +795,10 @@ def make_rupture_gdf(
     subfaults = fault_network['subfaults']
     rupture_df = fault_network[rup_df_key]
     sf_meshes = make_sf_rupture_meshes(
-        single_rup_df['patches'], single_rup_df['fault'], subfaults
+        single_rup_df['patches'],
+        single_rup_df['fault'],
+        subfaults,
+        same_size_arrays=same_size_arrays,
     )
     # converting to surfaces because get_boundary_3d doesn't take meshes
     sf_surfs = [SimpleFaultSurface(sf_mesh) for sf_mesh in sf_meshes]
@@ -806,7 +819,9 @@ def make_rupture_gdf(
     return rupture_gdf
 
 
-def merge_meshes_no_overlap(arrays, positions) -> np.ndarray:
+def merge_meshes_no_overlap(
+    arrays, positions, same_size_arrays: bool = True
+) -> np.ndarray:
     """
     Merges a list of arrays into a single array, with no overlap between
     the arrays.
@@ -825,9 +840,24 @@ def merge_meshes_no_overlap(arrays, positions) -> np.ndarray:
         Merged array.
     """
     # Check that all arrays have the same shape
-    first_shape = arrays[0].shape
-    for arr in arrays:
-        assert arr.shape == first_shape, "All arrays must have the same shape"
+    # Optional, but should be used for simple faults
+    # Kite faults can have different shapes
+    if same_size_arrays:
+        first_shape = arrays[0].shape
+        for arr in arrays:
+            assert (
+                arr.shape == first_shape
+            ), "All arrays must have the same shape"
+
+    else:
+        # check to see that all arrays share the same rows or same columns
+        row_lengths = [arr.shape[0] for arr in arrays]
+        col_lengths = [arr.shape[1] for arr in arrays]
+        assert (
+            len(set(row_lengths)) == 1 or len(set(col_lengths)) == 1
+        ), "All arrays must have the same number of rows or columns"
+        # `first_shape` is, in this case, the largest array
+        first_shape = (max(row_lengths), max(col_lengths))
 
     # check that all positions are unique and accounted for
     all_rows = sorted(list(set(pos[0] for pos in positions)))
@@ -861,7 +891,9 @@ def merge_meshes_no_overlap(arrays, positions) -> np.ndarray:
     return final_array
 
 
-def make_mesh_from_subfaults(subfaults: list[dict]) -> RectangularMesh:
+def make_mesh_from_subfaults(
+    subfaults: list[dict], same_size_arrays: bool = True
+) -> RectangularMesh:
     """
     Makes a RectangularMesh from a list of subfaults.
 
@@ -881,21 +913,26 @@ def make_mesh_from_subfaults(subfaults: list[dict]) -> RectangularMesh:
     big_lons = merge_meshes_no_overlap(
         [sf['surface'].mesh.lons for sf in subfaults],
         [sf['fault_position'] for sf in subfaults],
+        same_size_arrays=same_size_arrays,
     )
 
     big_lats = merge_meshes_no_overlap(
         [sf['surface'].mesh.lats for sf in subfaults],
         [sf['fault_position'] for sf in subfaults],
+        same_size_arrays=same_size_arrays,
     )
     big_depths = merge_meshes_no_overlap(
         [sf['surface'].mesh.depths for sf in subfaults],
         [sf['fault_position'] for sf in subfaults],
+        same_size_arrays=same_size_arrays,
     )
 
     return RectangularMesh(big_lons, big_lats, big_depths)
 
 
-def make_sf_rupture_mesh(rupture_indices, subfaults) -> RectangularMesh:
+def make_sf_rupture_mesh(
+    rupture_indices, subfaults, same_size_arrays: bool = True
+) -> RectangularMesh:
     """
     Makes a single-fault rupture mesh from a list of subfaults. This is
     a contiguous surface, unlike a multi-fault rupture surface.
@@ -913,12 +950,12 @@ def make_sf_rupture_mesh(rupture_indices, subfaults) -> RectangularMesh:
         Mesh composed of the meshes from all the subfaults in the rupture.
     """
     subs = [subfaults[i] for i in rupture_indices]
-    mesh = make_mesh_from_subfaults(subs)
+    mesh = make_mesh_from_subfaults(subs, same_size_arrays=same_size_arrays)
     return mesh
 
 
 def make_sf_rupture_meshes(
-    all_rupture_indices, faults, all_subfaults
+    all_rupture_indices, faults, all_subfaults, same_size_arrays: bool = True
 ) -> list[RectangularMesh]:
     """
     Makes a list of rupture meshes from a list of single-fault ruptures.
@@ -944,7 +981,9 @@ def make_sf_rupture_meshes(
     for i, rup_indices in enumerate(all_rupture_indices):
         try:
             subs_for_fault = grouped_subfaults[faults[i]]
-            mesh = make_sf_rupture_mesh(rup_indices, subs_for_fault)
+            mesh = make_sf_rupture_mesh(
+                rup_indices, subs_for_fault, same_size_arrays=same_size_arrays
+            )
             rup_meshes.append(mesh)
         except IndexError as e:
             logging.error(f"Problems with rupture {i}: " + str(e))
@@ -979,10 +1018,13 @@ def shapely_multipoly_to_geojson(multipoly, return_type='coords'):
 def export_ruptures_new(
     fault_network, rup_df_key='rupture_df_keep', outfile=None
 ):
-    subfault_gdf = make_subfault_gdf(fault_network['subfault_df'])
-    rupture_gdf = make_rupture_gdf(
-        fault_network, rup_df_key=rup_df_key, keep_sequences=True
-    )
+    # subfault_gdf = make_subfault_gdf(fault_network['subfault_df'])
+    if rup_df_key != 'rupture_gdf':
+        rupture_gdf = make_rupture_gdf(
+            fault_network, rup_df_key=rup_df_key, keep_sequences=True
+        )
+    else:
+        rupture_gdf = fault_network['rupture_gdf']
 
     outfile_type = outfile.split('.')[-1]
 
@@ -998,7 +1040,7 @@ def export_ruptures_new(
         features = []
         for i, rj in rup_json.items():
             f = geoms[i]
-            f["properties"] = {k: v for k, v in rj.items()}
+            f["properties"] = {k: v for k, v in rj.items() if k != 'geometry'}
             features.append(f)
 
         out_geojson = {"type": "FeatureCollection", "features": features}