datagenerator.util

   1import os
   2import numpy as np
   3
   4
   5def find_line_with_most_voxels(volume, voxel_thresholds, cfg):
   6    """
   7    Find the line with the most voxels
   8    ----------------------------------
   9
  10    Finds the line wtih the most voxels.
  11
  12    Loop over a volume and locate the inline which has the highest number of
  13    voxels which meet the threshold criteria
  14
  15    If multiple thresholds are used, first is used as a lower bound,
  16    and second is used as an upper bound
  17
  18    Parameters
  19    ----------
  20    volume : np.ndarray
  21        3D volume (ndarray)
  22    voxel_thresholds: list
  23        Used for thresholding
  24    cfg : dict
  25        model parameters
  26
  27    Returns
  28    -------
  29    inline_index : int
  30        Inline containing highest number of voxels meeting criteria.
  31    """
  32    max_voxels = 0
  33    inline_index = int(cfg.cube_shape[0] / 2)
  34    for ii in range(volume.shape[0]):
  35        voxels = volume[ii, ...].copy()
  36        try:  # if two values provided...
  37            voxels = np.where(
  38                ((voxels >= voxel_thresholds[0] & voxels < voxel_thresholds[1])),
  39                1.0,
  40                0.0,
  41            )
  42        except TypeError:  # Just use 1 threshold
  43            voxels = np.where(voxels > voxel_thresholds, 1.0, 0.0)
  44        if voxels[voxels == 1.0].size > max_voxels:
  45            max_voxels = voxels[voxels == 1.0].size
  46            inline_index = int(ii)
  47
  48    return inline_index
  49
  50
  51def plot_voxels_not_in_regular_layers(
  52    volume: np.ndarray,
  53    threshold: float,
  54    title: str,
  55    png_name: str,
  56    cfg
  57) -> None:
  58    """
  59        Plot voxels not in regular layers
  60
  61        Analyze voxel values not in regular layers by plotting a
  62        histogram of voxels above a given threshold.
  63
  64        Parameters
  65        ----------
  66        volume : np.ndarray
  67            3D volume (ndarray)
  68        threshold : float
  69            Threshold for voxel values.
  70        title : str
  71            Title for plot.
  72        png_name : str
  73            Name of png file to save.
  74        cfg : dict
  75            Model configurations.
  76
  77        Returns
  78        -------
  79        None
  80    """
  81    plt = import_matplotlib()
  82    voxels = volume[volume > threshold]
  83    plt.figure(1, figsize=(20, 15))
  84    plt.clf()
  85    plt.title(title)
  86    plt.hist(voxels, 81)
  87    plt.savefig(os.path.join(cfg.work_subfolder, png_name), format="png")
  88    plt.close(1)
  89
  90
  91def plot_xsection(volume, maps, line_num, title, png_name, cfg, cmap="prism") -> None:
  92    """
  93    Plot cross section.
  94
  95    Parameters
  96    ----------
  97    volume : np.ndarray
  98        The 3D volume as a numpy array
  99    maps : _type_
 100        _description_
 101    line_num : _type_
 102        _description_
 103    title : _type_
 104        _description_
 105    png_name : _type_
 106        _description_
 107    cfg : _type_
 108        _description_
 109    cmap : str, optional
 110        _description_, by default "prism"
 111    
 112    Returns
 113    -------
 114    None
 115    """
 116    plt = import_matplotlib()
 117    plt.clf()
 118    plt.title(f"{title}\nInline: {line_num}")
 119    plt.figure(1, figsize=(20, 15))
 120    # If depth maps is infilled and volume is not, update the Z axis of the volume
 121    # Pick the 75%'th horizon to check
 122    if (
 123        np.max(maps[:, :, -int(maps.shape[-1] / 4)]) > volume.shape[-1]
 124    ):  # don't use last horizon
 125        plt.imshow(
 126            np.fliplr(
 127                np.rot90(
 128                    volume[line_num, :, : (volume.shape[-1] * cfg.infill_factor) - 1], 3
 129                )
 130            ),
 131            aspect="auto",
 132            cmap=cmap,
 133        )
 134    else:
 135        plt.imshow(
 136            np.fliplr(np.rot90(volume[line_num, ...], 3)), aspect="auto", cmap=cmap
 137        )
 138    plt.colorbar()
 139    plt.ylim((volume.shape[-1], 0))
 140    for i in range(0, maps.shape[-1], 1):
 141        plt.plot(range(cfg.cube_shape[0]), maps[line_num, :, i], "k-", lw=0.3)
 142    plt.savefig(os.path.join(cfg.work_subfolder, png_name), format="png")
 143    plt.close(1)
 144
 145
 146def plot_3D_faults_plot(cfg, faults, plot_faults=True, plot_throws=True) -> None:
 147    """
 148    Plot 3D faults plot
 149    -------------------
 150
 151    Plots the faults in a 3D plot.
 152
 153    Parameters
 154    ----------
 155    cfg : _type_
 156        The configuration.
 157    faults : np.ndarray
 158        The faults a numpy object.
 159    plot_faults : bool, optional
 160        Whether to plot the faults or not, by default True
 161    plot_throws : bool, optional
 162        Whether to plot the fault throws or not, by default True
 163    
 164    Returns
 165    -------
 166    None
 167    """
 168    from plotly.offline import plot
 169    import plotly.graph_objects as go
 170
 171    fi1 = faults.fault_planes
 172    decimation_factor = 2
 173    x1, y1, z1 = np.where(
 174        fi1[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5
 175    )
 176    x1 *= decimation_factor
 177    y1 *= decimation_factor
 178    z1 *= -decimation_factor
 179    trace1 = go.Scatter3d(
 180        x=x1,
 181        y=y1,
 182        z=z1,
 183        name="fault_segments",
 184        mode="markers",
 185        marker=dict(size=2.0, color="blue", opacity=0.025),
 186    )
 187
 188    if plot_faults:
 189        fi2 = faults.fault_intersections
 190        x2, y2, z2 = np.where(
 191            fi2[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.25
 192        )
 193        x2 *= decimation_factor
 194        y2 *= decimation_factor
 195        z2 *= -decimation_factor
 196        trace2 = go.Scatter3d(
 197            x=x2,
 198            y=y2,
 199            z=z2,
 200            name="fault_intersection_segments",
 201            mode="markers",
 202            marker=dict(
 203                size=1.5,
 204                color="red",  # set color to an array/list of desired values
 205                colorscale="Viridis",  # choose a colorscale
 206                opacity=0.025,
 207            ),
 208        )
 209
 210        data = [trace1, trace2]
 211        layout = go.Layout(
 212            title=go.layout.Title(
 213                text="Fault Segments & Intersections"
 214                + "<br>"
 215                + '<span style="font-size: 12px;">'
 216                + cfg.work_subfolder
 217                + "</span>",
 218                xref="paper",
 219                x=0,
 220            )
 221        )
 222        camera = dict(eye=dict(x=1.25 / 1.2, y=1.25 / 1.2, z=0.75 / 1.2))
 223        fig3 = go.Figure(data=data, layout=layout)
 224        fig3.update_layout(
 225            scene_camera=camera,
 226            autosize=False,
 227            width=960,
 228            height=720,
 229            scene=dict(
 230                aspectmode="manual",
 231                aspectratio=dict(x=1, y=1, z=1),
 232                xaxis=dict(
 233                    nticks=4,
 234                    range=[0, 300],
 235                ),
 236                yaxis=dict(
 237                    nticks=4,
 238                    range=[0, 300],
 239                ),
 240                zaxis=dict(
 241                    nticks=4,
 242                    range=[-1260, 0],
 243                ),
 244            ),
 245            margin=dict(r=20, l=10, b=10, t=25),
 246        )
 247        output_filename = os.path.join(
 248            cfg.work_subfolder,
 249            f"qcplot_fault_segments_and_intersections_3D__{cfg.date_stamp}",
 250        )
 251        plot(fig3, filename=f"{output_filename}.html", auto_open=False)
 252        try:
 253            fig3.write_image(f"{output_filename}.png")
 254        except:
 255            print("png file saving failed")
 256        print(
 257            "\n   ... util/plot_3D_faults_plot finished creating 3D plot at:\n       "
 258            + output_filename
 259        )
 260
 261    if plot_throws:
 262        fi3 = faults.fault_plane_throw
 263        z3 = fi3[x1, y1, -z1]
 264        trace3 = go.Scatter3d(
 265            x=x1,
 266            y=y1,
 267            z=z1,
 268            name="fault_intersection_segments",
 269            mode="markers",
 270            marker=dict(
 271                size=4.0,
 272                color=z3,  # set color to an array/list of desired values
 273                colorscale="Viridis",  # choose a colorscale
 274                opacity=0.05,
 275            ),
 276        )
 277        data = [trace3]
 278        layout = go.Layout(
 279            title=go.layout.Title(
 280                text="Fault Throw along fault plane"
 281                + "<br>"
 282                + '<span style="font-size: 12px;">'
 283                + cfg.work_subfolder
 284                + "</span>",
 285                xref="paper",
 286                x=0,
 287            )
 288        )
 289        camera = dict(eye=dict(x=1.25 / 1.2, y=1.25 / 1.2, z=0.75 / 1.2))
 290        fig4 = go.Figure(data=data, layout=layout)
 291        fig4.update_layout(
 292            scene_camera=camera,
 293            autosize=False,
 294            width=960,
 295            height=720,
 296            scene=dict(
 297                aspectmode="manual",
 298                aspectratio=dict(x=1, y=1, z=1),
 299                xaxis=dict(
 300                    nticks=4,
 301                    range=[0, 300],
 302                ),
 303                yaxis=dict(
 304                    nticks=4,
 305                    range=[0, 300],
 306                ),
 307                zaxis=dict(
 308                    nticks=4,
 309                    range=[-1260, 0],
 310                ),
 311            ),
 312            margin=dict(r=20, l=10, b=10, t=25),
 313        )
 314        output_filename = os.path.join(
 315            cfg.work_subfolder, f"qcplot_fault_segments_throw_3D__{cfg.date_stamp}"
 316        )
 317        plot(fig4, filename=f"{output_filename}.html", auto_open=False)
 318        try:
 319            fig4.write_image(f"{output_filename}.png")
 320        except:
 321            print("png file saving failed")
 322        print(
 323            "\n   ... util/plot_3D_faults_plot finished creating 3D plot at:\n       "
 324            + output_filename
 325        )
 326
 327
 328def plot_3D_closure_plot(cfg, closures, plot_closures=True, plot_strat_closures=True) -> None:
 329    """
 330    Plot 3D closures
 331    ----------------
 332
 333    Plot the closures in 3D.
 334
 335    Parameters
 336    ----------
 337    cfg : dict
 338        The configuration.
 339    closures : np.ndarray
 340        Closures array.
 341    plot_closures : bool, optional
 342        Whether to plot the closures or not, by default True
 343    plot_strat_closures : bool, optional
 344        Whether to plot statigraphic closures or not, by default True
 345    
 346    Returns
 347    -------
 348    None
 349    """
 350    from plotly.offline import plot
 351    import plotly.graph_objects as go
 352
 353    fi1 = closures.gas_closures
 354    decimation_factor = 2
 355    x1, y1, z1 = np.where(
 356        fi1[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5
 357    )
 358    x1 *= decimation_factor
 359    y1 *= decimation_factor
 360    z1 *= -decimation_factor
 361
 362    fi2 = closures.oil_closures
 363    x2, y2, z2 = np.where(
 364        fi2[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.25
 365    )
 366    x2 *= decimation_factor
 367    y2 *= decimation_factor
 368    z2 *= -decimation_factor
 369
 370    fi3 = closures.brine_closures
 371    x3, y3, z3 = np.where(
 372        fi3[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5
 373    )
 374    x3 *= decimation_factor
 375    y3 *= decimation_factor
 376    z3 *= -decimation_factor
 377
 378    size = 2.0
 379    opacity = 0.15
 380
 381    if plot_closures:
 382        trace1 = go.Scatter3d(
 383            x=x1,
 384            y=y1,
 385            z=z1,
 386            name="gas_segments",
 387            mode="markers",
 388            marker=dict(size=size, color="red", opacity=opacity),
 389        )
 390
 391        trace2 = go.Scatter3d(
 392            x=x2,
 393            y=y2,
 394            z=z2,
 395            name="oil_segments",
 396            mode="markers",
 397            marker=dict(
 398                size=size,
 399                color="green",  # set color to an array/list of desired values
 400                colorscale="Viridis",  # choose a colorscale
 401                opacity=opacity,
 402            ),
 403        )
 404
 405        trace3 = go.Scatter3d(
 406            x=x3,
 407            y=y3,
 408            z=z3,
 409            name="brine_segments",
 410            mode="markers",
 411            marker=dict(
 412                size=size,
 413                color="blue",  # set color to an array/list of desired values
 414                colorscale="Viridis",  # choose a colorscale
 415                opacity=opacity / 5.0,
 416            ),
 417        )
 418
 419        data = [trace1, trace2, trace3]
 420        layout = go.Layout(
 421            title=go.layout.Title(
 422                text="Closure Segments (Gas, Oil, Brine)"
 423                + "<br>"
 424                + '<span style="font-size: 12px;">'
 425                + cfg.work_subfolder
 426                + "</span>",
 427                xref="paper",
 428                x=0,
 429            )
 430        )
 431        camera = dict(eye=dict(x=1.25 / 1.2, y=1.25 / 1.2, z=0.75 / 1.2))
 432        fig4 = go.Figure(data=data, layout=layout)
 433        fig4.update_layout(
 434            scene_camera=camera,
 435            autosize=False,
 436            width=960,
 437            height=720,
 438            scene=dict(
 439                aspectmode="manual",
 440                aspectratio=dict(x=1, y=1, z=1),
 441                xaxis=dict(
 442                    nticks=4,
 443                    range=[0, 300],
 444                ),
 445                yaxis=dict(
 446                    nticks=4,
 447                    range=[0, 300],
 448                ),
 449                zaxis=dict(
 450                    nticks=4,
 451                    range=[-1260, 0],
 452                ),
 453            ),
 454            # width=650,
 455            margin=dict(r=20, l=10, b=10, t=25),
 456        )
 457        output_filename = os.path.join(
 458            cfg.work_subfolder, f"qcplot_closure_segments_3D_{cfg.date_stamp}"
 459        )
 460        plot(fig4, filename=f"{output_filename}.html", auto_open=False)
 461        try:
 462            fig4.write_image(f"{output_filename}.png")
 463        except ValueError:
 464            print("png file saving failed")
 465        print(
 466            "\n   ... util/plot_3D_closure_plot finished creating 3D plot at:\n       "
 467            + output_filename
 468        )
 469
 470    if plot_strat_closures:
 471        ###-------------------------------------------------------------------------
 472        ### Make 3D plot of strat traps (onlap/pinch-out traps)
 473        ###-------------------------------------------------------------------------
 474        fi4 = closures.strat_closures
 475        decimation_factor = 2
 476        x4, y4, z4 = np.where(
 477            (fi4[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 478            & (fi1[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 479        )
 480        x4 *= decimation_factor
 481        y4 *= decimation_factor
 482        z4 *= -decimation_factor
 483
 484        x5, y5, z5 = np.where(
 485            (fi4[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 486            & (fi2[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 487        )
 488        x5 *= decimation_factor
 489        y5 *= decimation_factor
 490        z5 *= -decimation_factor
 491
 492        x6, y6, z6 = np.where(
 493            (fi4[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 494            & (fi3[::decimation_factor, ::decimation_factor, ::decimation_factor] > 0.5)
 495        )
 496        x6 *= decimation_factor
 497        y6 *= decimation_factor
 498        z6 *= -decimation_factor
 499
 500        size = 2.0
 501        opacity = 0.15
 502
 503        trace1 = go.Scatter3d(
 504            x=x4,
 505            y=y4,
 506            z=z4,
 507            name="gas_strat",
 508            mode="markers",
 509            marker=dict(size=size, color="red", opacity=opacity),
 510        )
 511
 512        trace2 = go.Scatter3d(
 513            x=x5,
 514            y=y5,
 515            z=z5,
 516            name="oil_strat",
 517            mode="markers",
 518            marker=dict(
 519                size=size,
 520                color="green",  # set color to an array/list of desired values
 521                colorscale="Viridis",  # choose a colorscale
 522                opacity=opacity,
 523            ),
 524        )
 525
 526        trace3 = go.Scatter3d(
 527            x=x6,
 528            y=y6,
 529            z=z6,
 530            name="brine_strat",
 531            mode="markers",
 532            marker=dict(
 533                size=size,
 534                color="blue",  # set color to an array/list of desired values
 535                colorscale="Viridis",  # choose a colorscale
 536                opacity=opacity / 5.0,
 537            ),
 538        )
 539
 540        data = [trace1, trace2, trace3]
 541        layout = go.Layout(
 542            title=go.layout.Title(
 543                text="Stratigraphic Closure Segments (Gas, Oil, Brine)"
 544                + "<br>"
 545                + '<span style="font-size: 12px;">'
 546                + cfg.work_subfolder
 547                + "</span>",
 548                xref="paper",
 549                x=0,
 550            )
 551        )
 552        camera = dict(eye=dict(x=1.25 / 1.2, y=1.25 / 1.2, z=0.75 / 1.2))
 553        fig5 = go.Figure(data=data, layout=layout)
 554        fig5.update_layout(
 555            scene_camera=camera,
 556            autosize=False,
 557            width=960,
 558            height=720,
 559            scene=dict(
 560                aspectmode="manual",
 561                aspectratio=dict(x=1, y=1, z=1),
 562                xaxis=dict(
 563                    nticks=4,
 564                    range=[0, 300],
 565                ),
 566                yaxis=dict(
 567                    nticks=4,
 568                    range=[0, 300],
 569                ),
 570                zaxis=dict(
 571                    nticks=4,
 572                    range=[-1260, 0],
 573                ),
 574            ),
 575            # width=650,
 576            margin=dict(r=20, l=10, b=10, t=25),
 577        )
 578        output_filename = os.path.join(
 579            cfg.work_subfolder, f"qcplot_strat_closure_segments_3D_{cfg.date_stamp}"
 580        )
 581        plot(fig5, filename=f"{output_filename}.html", auto_open=False)
 582        try:
 583            fig5.write_image(f"{output_filename}.png")
 584        except:
 585            print("png file saving failed")
 586        print(
 587            "\n   ... util/plot_3D_strat_closure_plot finished creating 3D plot at:\n       "
 588            + output_filename
 589        )
 590
 591
 592def plot_facies_qc(cfg, lith, seismic, facies, maps) -> None:
 593    """
 594    Plot Facies for QC
 595    ------------------
 596
 597    Plots the facies to aid in QC of the model.
 598
 599    Parameters
 600    ----------
 601    cfg : dict
 602        The configuration used in the model.
 603    lith : np.ndarray
 604        The lithology model.
 605    seismic : np.ndarray
 606        The seismic model.
 607    facies : np.ndarray
 608        The facies model.
 609    maps : dict
 610        The maps used in the model.
 611    
 612    Returns
 613    -------
 614    None
 615    """
 616    plt = import_matplotlib()
 617    from itertools import groupby
 618    import matplotlib as mpl
 619
 620    # Plot Centre Inline of facies and seismic
 621    iline = lith.shape[0] // 2
 622    avg_sand_unit_thickness = np.mean(
 623        [b for a, b in [(k, sum(1 for i in g)) for k, g in groupby(facies)] if a == 1.0]
 624    )
 625    textstr = "\n".join(
 626        (
 627            f"Sand % Input: {cfg.sand_layer_pct:.1%}",
 628            f"Sand % Actual: {facies[facies == 1.].size / facies.size:.1%}",
 629            f"Sand thickness (layers) Input: {cfg.sand_layer_thickness}",
 630            f"Sand thickness (layers) Actual: {avg_sand_unit_thickness:.1f}",
 631            f"Sand Voxel % in model: {lith[lith[:] == 1].size / lith[lith[:] >= 0].size:.1%}",
 632        )
 633    )
 634    fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(20, 15), sharey=True)
 635    fig.suptitle("Facies QC Plot")
 636    axs[0].set_ylim(cfg.cube_shape[-1])
 637    if cfg.include_salt and np.max(lith[iline, ...]) > 1:
 638        axs[0].imshow(
 639            np.fliplr(np.rot90(lith[iline, ...], 3)),
 640            aspect="auto",
 641            cmap=mpl.colors.ListedColormap(["blue", "saddlebrown", "gold", "grey"]),
 642        )
 643    else:
 644        axs[0].imshow(
 645            np.fliplr(np.rot90(lith[iline, ...], 3)),
 646            aspect="auto",
 647            cmap=mpl.colors.ListedColormap(["blue", "saddlebrown", "gold"]),
 648        )
 649    _img = axs[1].imshow(
 650        np.fliplr(np.rot90(seismic[iline, ...], 3)),
 651        aspect="auto",
 652        cmap="Greys",
 653        vmin=-300,
 654        vmax=300,
 655    )
 656    # fig.colorbar(_img, ax=axs[1])
 657    props = dict(boxstyle="round", alpha=0.5)
 658    # Add textbox with textstr to facies subplot
 659    axs[0].text(
 660        0.05,
 661        0.95,
 662        textstr,
 663        transform=axs[0].transAxes,
 664        fontsize=14,
 665        verticalalignment="top",
 666        bbox=props,
 667    )
 668    for i in range(maps.shape[-1]):
 669        axs[0].plot(range(cfg.cube_shape[0]), maps[iline, :, i], "k-", lw=0.3)
 670
 671    fig.savefig(os.path.join(cfg.work_subfolder, "QC_plot__Facies_FullStackCumulativeSum.png"))
 672    plt.close()
 673
 674
 675def infill_surface(surface: np.ndarray) -> np.ndarray:
 676    """
 677    Infill Surface
 678    --------------
 679
 680    Infill surfaces.
 681
 682    Fill holes in input 2D surface.
 683    Holes have value of either nan or zero.
 684    Fill using replacement with interpolated (nearest-neighbor) value.
 685
 686    Parameters
 687    ----------
 688    surface : np.ndarray
 689        2D numpy array
 690
 691    Returns
 692    -------
 693    surface_infilled : np.ndarray
 694        (2D) surface_infilled
 695    """
 696    from scipy.interpolate import NearestNDInterpolator
 697
 698    # get x,y,z indices for non-nan points on surface
 699    xx, yy = np.meshgrid(
 700        range(surface.shape[0]), range(surface.shape[1]), sparse=False, indexing="ij"
 701    )
 702    x_no_nan = xx[~np.isnan(surface)]
 703    y_no_nan = yy[~np.isnan(surface)]
 704    z_surface = surface[~np.isnan(surface)]
 705
 706    if z_surface.flatten().shape[0] > 2:
 707        # set up nearest-neighbor interpolator function
 708        xy = zip(x_no_nan, y_no_nan)
 709        nn = NearestNDInterpolator(xy, z_surface)
 710        # interpolate at every x,y on regular grid
 711        surface_nn = nn(xx, yy)
 712        # replace input surface with interpolated (nearest-neighbor)
 713        # - if input value is either nan or zero
 714        surface_infilled = surface.copy()
 715        surface_infilled[np.isnan(surface)] = surface_nn[np.isnan(surface)]
 716        surface_infilled[surface == 0] = surface_nn[surface == 0]
 717    else:
 718        surface_infilled = z_surface.copy()
 719
 720    if surface_infilled[np.isnan(surface_infilled)].shape[0] > 0:
 721        count = surface_infilled[np.isnan(surface_infilled)].shape[0]
 722        print(
 723            f"\n\n\n   ...inside infill_surface: there are some NaN values in the surface. count = {count}"
 724        )
 725
 726    return surface_infilled
 727
 728
 729def mute_above_seafloor(surface, xyz):
 730    """
 731    Mute data above seafloor
 732    ------------------------
 733
 734    Mute a cube above a surface that contains
 735    indices for the 3rd axis of the cube.
 736
 737    Parameters
 738    ----------
 739    surface : np.ndarray
 740        Mute above this (2D) surface
 741    xyz : np.ndarray
 742        Cube to which muting is applied.
 743
 744    Returns
 745    -------
 746    xyz_muted : np.ndarray
 747        Muted (3D array)
 748    """
 749
 750    krange = np.arange(xyz.shape[2])
 751
 752    # broadcast indices of vertical axis to 3D
 753    krange = krange.reshape((1, 1, len(krange)))
 754
 755    if np.isnan(np.min(surface)) or 0 in surface and not np.all(surface == 0):
 756        # If surface contains nan's or at least one 0, this will return True, therefore...
 757        # fill empty picks with nearest-neighbor infilling
 758        # If all 0's, skip it
 759        surface_infilled = infill_surface(surface)
 760    else:
 761        surface_infilled = surface
 762
 763    # broadcast surface to 3D
 764    surface_3d = surface_infilled.reshape((surface.shape[0], surface.shape[1], 1))
 765
 766    # apply muting to copy of xyz
 767    xyz_muted = xyz.copy()
 768    xyz_muted[krange < surface_3d] *= 0.0
 769
 770    return xyz_muted
 771
 772
 773def is_it_in_hull(hull, p) -> np.ndarray[bool]:
 774    """
 775    Is it in hull?
 776    --------------
 777    Checks if point is in hull
 778
 779    Test if points in `p` are in `hull`
 780
 781    `p` should be a `NxK` coordinates of `N` points in `K` dimensions
 782    `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the
 783    coordinates of `M` points in `K`dimensions for which Delaunay triangulation
 784    will be computed
 785
 786    Parameters
 787    ----------
 788    hull : np.ndarray
 789        The hull object.
 790    p : np.ndarray
 791        Point(s) to check.
 792
 793    Returns
 794    -------
 795    np.ndarray[bool]
 796        An array containing the booleans where the object is hull.
 797
 798
 799    """
 800    from scipy.spatial import Delaunay
 801
 802    if not isinstance(hull, Delaunay):
 803        hull = Delaunay(hull)
 804
 805    return hull.find_simplex(p) >= 0
 806
 807
 808# Utility Functions
 809def even_odd(number):
 810    if number % 2 == 0:
 811        return "Even"
 812    else:
 813        return "Odd"
 814
 815
 816def next_odd(number):
 817    if even_odd(number) == "Even":
 818        odd_number = number + 1
 819    else:
 820        odd_number = number
 821    return odd_number
 822
 823
 824def find_index_next_bigger(mylist, mynumber):
 825    # return the index of the next number in mylist that is bigger (or equal to) mynumber
 826    # example:   mynumber=46;  mylist=[1,5,13,25,36,48,55,67,73,85,90]
 827    # example:  find_index_next_bigger(mylist,mynumber) returns 5
 828    return next((mylist.index(n) for n in mylist if n >= mynumber), len(mylist))
 829
 830
 831def linearinfill(x, y, newx):
 832    # ***************************************************************************************
 833    #
 834    #   Function to return data spaced at locations specified by input variable 'newx' after
 835    #   fitting linear function. Input data specified by 'x' and 'y' so data
 836    #   can be irregularly sampled.
 837    #
 838    # ***************************************************************************************
 839    s = np.interp(newx, x, y)
 840    return s
 841
 842
 843def check_resident_memory():
 844    import os
 845    import sys
 846
 847    if "linux" in sys.platform:
 848        _proc_status = "/proc/%d/status" % os.getpid()
 849        _scale = {
 850            "kB": 1024.0,
 851            "mB": 1024.0 * 1024.0,
 852            "KB": 1024.0,
 853            "MB": 1024.0 * 1024.0,
 854        }
 855
 856        t = open(_proc_status)
 857        v = t.read()
 858        t.close()
 859
 860        i = v.index("VmRSS:")
 861        v = v[i:].split(None, 3)
 862
 863        return (float(v[1]) * _scale[v[2]]) / 1024 ** 3
 864    else:
 865        return 0.0
 866
 867
 868def import_matplotlib():
 869    """
 870    Function to perform matplotlib imports with non-interactive backend when running in background
 871    DPG 17/2/19
 872    """
 873    import sys
 874
 875    if bool(getattr(sys, "ps1", sys.flags.interactive)):
 876        # this is an interactive session
 877        from matplotlib import pyplot as plt
 878
 879        plt.ion()
 880    else:
 881        # import matplotlib with non-interactive backend ('Agg')
 882        import matplotlib
 883
 884        matplotlib.use("Agg")
 885        from matplotlib import pyplot as plt
 886
 887        plt.ioff()
 888    return plt
 889
 890
 891def write_data_to_hdf(n, d, h5file):
 892    import h5py
 893
 894    with h5py.File(h5file, "a") as hf:
 895        hf.create_dataset(name=n, data=d, compression="lzf")
 896
 897
 898# Useful functions, not directly used for data generation
 899def qc_folders(directory):
 900    """Count how many model_parameter.txt files contain elapse time
 901
 902    Use this to remove models that failed to complete for whatever reason
 903    """
 904    import glob
 905    import os
 906
 907    valid_models = list()
 908    model_closures = dict()
 909    # Grep for elapse in model_parameters.txt files in geocrawler subdirectories
 910    par = glob.glob(os.path.join(directory, "geoc*/model_parameters*.txt"))
 911    # Elapse time
 912    for p in par:
 913        with open(p, "r") as f:
 914            for line in f.readlines():
 915                if "elapse" in line:
 916                    valid_models.append(os.path.dirname(p))
 917                if "Number of Closures" in line:
 918                    num_closures = line.split()[-1]
 919                    model_closures[os.path.dirname(p)] = int(num_closures)
 920
 921    # Invalid folders
 922    invalid_models = list(
 923        set(glob.glob("{}/geoc*".format(directory))) - set(valid_models)
 924    )
 925    # Remove tar.gz files from invalid_models
 926    invalid_final = [x for x in invalid_models if not x.endswith("tar.gz")]
 927
 928    # Also sort closures by number of closures
 929    sorted_closures = [
 930        (k, model_closures[k])
 931        for k in sorted(model_closures, key=model_closures.get, reverse=True)
 932    ]
 933
 934    return invalid_final, valid_models, sorted_closures
 935
 936
 937def hanflat(inarray, pctflat):
 938
 939    # ***************************************************************************************
 940    #
 941    #   Function applies a Hanning taper to ends of "inarray".
 942    #   Center "pctflat" of samples remain unchanged.
 943    #
 944    #   Parameters:
 945    #   array :       array of values to have ends tapered
 946    #   pctflat :     this percent of  samples to remain unchanged (e.g. 0.80)
 947    #
 948    #   Returns weights (not weights applied to input array)
 949    #
 950    # ***************************************************************************************
 951
 952    numsamples = len(inarray)
 953    lowflatindex = int(round(numsamples * (1.0 - pctflat) / 2.0))
 954    hiflatindex = numsamples - lowflatindex
 955
 956    # get hanning for numsamples*(1.0-pctflat)
 957    hanwgt = np.hanning(len(inarray) - (hiflatindex - lowflatindex))
 958
 959    # piece together hanning weights at ends and weights=1.0 in flat portion
 960    outarray = np.ones(len(inarray), dtype=float)
 961    outarray[:lowflatindex] = hanwgt[:lowflatindex]
 962    outarray[hiflatindex:] = hanwgt[numsamples - hiflatindex :]
 963
 964    return outarray
 965
 966
 967def hanflat2d(inarray2d, pctflat):
 968    hanwgt1 = hanflat(np.ones((inarray2d.shape[0])), pctflat)
 969    hanwgt2 = hanflat(np.ones((inarray2d.shape[1])), pctflat)
 970    return np.sqrt(np.outer(hanwgt1, hanwgt2))
 971
 972
 973def hanflat3d(inarray3d, pctflat):
 974    hanwgt1 = hanflat(np.ones((inarray3d.shape[0])), pctflat)
 975    hanwgt2 = hanflat(np.ones((inarray3d.shape[1])), pctflat)
 976    hanwgt3 = hanflat(np.ones((inarray3d.shape[2])), pctflat)
 977    _a = np.outer(hanwgt1, hanwgt2)
 978    _b = np.multiply.outer(_a.ravel(), hanwgt3.ravel())
 979    return _b.reshape(inarray3d.shape) ** (1.0 / 3.0)
 980
 981
 982def create_3D_qc_plot(
 983    cfg,
 984    seismic,
 985    brine,
 986    oil,
 987    gas,
 988    fault_segments,
 989    title_text,
 990    camera_steps=20
 991) -> None:
 992    """
 993    Creates a 3D QC plot
 994    --------------------
 995
 996    Generates a 3D QC plot for aiding in QC.
 997
 998    Parameters
 999    ----------
1000    cfg : _type_
1001        _description_
1002    seismic : _type_
1003        _description_
1004    brine : _type_
1005        _description_
1006    oil : _type_
1007        _description_
1008    gas : _type_
1009        _description_
1010    fault_segments : _type_
1011        _description_
1012    title_text : _type_
1013        _description_
1014    camera_steps : int, optional
1015        _description_, by default 20
1016
1017    Returns
1018    -------
1019    None
1020    """
1021
1022    from scipy.ndimage import gaussian_filter1d, gaussian_filter
1023    from vtkplotter import addons, show, Text2D, Volume, Plotter
1024    from vtkplotter.pyplot import cornerHistogram
1025    from vtkplotter.vtkio import screenshot
1026
1027    # custom lighting for meshes
1028    ambient = 0.4
1029    diffuse = 0.1
1030    specular = 0.5
1031    specular_power = 30
1032
1033    # Brine
1034    # rescale the vertical axis to have similar height to x & y axes width
1035    v_exagg = 1.0 / (cfg.cube_shape[-1] // cfg.cube_shape[0])
1036    vol_br = Volume(brine[:, :, ::-1], spacing=(1, 1, v_exagg))
1037    vol_br = Volume(vol_br.GetMapper().GetInput())
1038    print(
1039        f" ... brine_closure amp range = {brine.min()}, {brine.mean()}, {brine.max()}"
1040    )
1041    brine_surface = (
1042        vol_br.isosurface(threshold=0.99 * brine.max())
1043        .c("blue")
1044        .lighting(
1045            ambient=ambient,
1046            diffuse=diffuse,
1047            specular=specular + 0.2,
1048            specularPower=specular_power,
1049        )
1050    )
1051    # Oil
1052    vol_oil = Volume(oil[:, :, ::-1], spacing=(1, 1, v_exagg))
1053    vol_oil = Volume(vol_oil.GetMapper().GetInput())
1054    print(f" ... oil_closure amp range = {oil.min()}, {oil.mean()}, {oil.max()}")
1055    oil_surface = (
1056        vol_oil.isosurface(threshold=0.99 * brine.max())
1057        .c("green")
1058        .lighting(
1059            ambient=ambient,
1060            diffuse=diffuse,
1061            specular=specular,
1062            specularPower=specular_power,
1063        )
1064    )
1065    # Gas
1066    vol_gas = Volume(gas[:, :, ::-1], spacing=(1, 1, v_exagg))
1067    vol_gas = Volume(vol_gas.GetMapper().GetInput())
1068    print(f" ... gas amp range = {gas.min()}, {gas.mean()}, {gas.max()}")
1069    gas_surface = (
1070        vol_gas.isosurface(threshold=0.99 * brine.max())
1071        .c("red")
1072        .lighting(
1073            ambient=ambient,
1074            diffuse=diffuse,
1075            specular=specular,
1076            specularPower=specular_power,
1077        )
1078    )
1079    # Faults
1080    try:
1081        sigma = 0.35
1082        threshold = 0.0759
1083        fault_segments = gaussian_filter(fault_segments, sigma=(sigma, sigma, 1))
1084        fault_segments = gaussian_filter1d(fault_segments, sigma=sigma, axis=2)
1085        vol_fault = Volume(fault_segments[:, :, ::-1], spacing=(1, 1, v_exagg))
1086        vol_fault = Volume(vol_fault.GetMapper().GetInput())
1087        fault_surface = (
1088            vol_fault.isosurface(threshold=threshold)
1089            .c((0.5, 0.5, 0.5))
1090            .alpha(0.05)
1091            .lighting(
1092                ambient=ambient,
1093                diffuse=diffuse,
1094                specular=specular * 1.5,
1095                specularPower=specular_power * 2,
1096            )
1097        )
1098        include_faults = True
1099    except:
1100        include_faults = False
1101
1102    # Seismic
1103    # re-order z so that big indices are deep
1104    vol = Volume(seismic[:, :, ::-1], spacing=(1, 1, v_exagg))
1105    vol.addScalarBar3D()
1106    vol = Volume(vol.GetMapper().GetInput())
1107
1108    # Create the 3D scene and save images
1109    la, ld = 0.9, 0.9  # ambient, diffuse
1110    cmaps = ["bone_r", "gist_ncar_r", "jet", "Spectral_r", "hot_r", "gist_earth_r"]
1111
1112    box = vol.box().wireframe().alpha(0)  # make an invisible box
1113
1114    vp = Plotter(
1115        axes=1,
1116        bg="white",
1117        bg2="lightblue",
1118        size=(900, 600),
1119        title=title_text,
1120        interactive=False,
1121        offscreen=True,
1122    )
1123    vp.show(box, newPlotter=True, interactive=False)
1124    # vp.showInset(vol, c=cmaps[0], pos=(1, 1), size=0.3, draggable=False)
1125
1126    # inits
1127    visibles = [None, None, None]
1128    cmap = cmaps[0]
1129    dims = vol.dimensions()
1130    i_init = 0
1131    msh = vol.xSlice(i_init).pointColors(cmap=cmap).lighting("", la, ld, 0)
1132    msh.addScalarBar(pos=(0.04, 0.0), horizontal=True, titleFontSize=0)
1133    vp.renderer.AddActor(msh)
1134    visibles[0] = msh
1135    j_init = 0
1136    msh2 = vol.ySlice(j_init).pointColors(cmap=cmap).lighting("", la, ld, 0)
1137    msh2.addScalarBar(pos=(0.04, 0.0), horizontal=True, titleFontSize=0)
1138    vp.renderer.AddActor(msh)
1139    visibles[1] = msh2
1140    k_init = 0
1141    msh3 = vol.zSlice(k_init).pointColors(cmap=cmap).lighting("", la, ld, 0, 0)
1142    msh3.addScalarBar(pos=(0.04, 0.0), horizontal=True, titleFontSize=0)
1143    vp.renderer.AddActor(msh3)
1144    visibles[2] = msh3
1145
1146    # the colormap button
1147    def buttonfunc():
1148        global cmap
1149        bu.switch()
1150        cmap = bu.status()
1151        for mesh in visibles:
1152            if mesh:
1153                mesh.pointColors(cmap=cmap)
1154        vp.renderer.RemoveActor(mesh.scalarbar)
1155        mesh.scalarbar = addons.addScalarBar(
1156            mesh, pos=(0.04, 0.0), horizontal=True, titleFontSize=0
1157        )
1158        vp.renderer.AddActor(mesh.scalarbar)
1159
1160    bu = vp.addButton(
1161        buttonfunc,
1162        pos=(0.27, 0.005),
1163        states=cmaps,
1164        c=["db"] * len(cmaps),
1165        bc=["lb"] * len(cmaps),
1166        size=14,
1167        bold=True,
1168    )
1169
1170    hist = cornerHistogram(
1171        seismic,
1172        s=0.2,
1173        bins=75,
1174        logscale=0,
1175        pos=(0.02, 0.02),
1176        c=((0.0, 0.0, 0.75)),
1177        bg=((0.1, 0.1, 0.1)),
1178        alpha=0.7,
1179    )
1180
1181    # Show the meshes and start interacting
1182    viewpoint1 = [-200.0, -200.0, 500.0]
1183    focalpoint1 = [300, 300, 150]
1184    vp.addLight(pos=viewpoint1, focalPoint=focalpoint1)
1185
1186    viewpoint2 = [-500.0, -500.0, 500.0]
1187    focalpoint2 = [150, 150, 150]
1188    vp.addLight(pos=viewpoint2, focalPoint=focalpoint2)
1189
1190    viewpoint3 = [-400.0, -400.0, 600.0]
1191    focalpoint3 = [150, 150, 150]
1192    vp.addLight(pos=viewpoint3, focalPoint=focalpoint3)
1193
1194    # Can now add any other object to the Plotter scene:
1195    text0 = Text2D(title_text, s=1.0, font="arial")
1196
1197    # set camera position
1198    start_pos = [825.0, 350.0, 600.0]
1199    start_vu = [-0.4, -0.1, 0.92]
1200    vp.camera.SetPosition(start_pos)
1201    vp.camera.SetFocalPoint([150.0, 150.0, 150.0])
1202    vp.camera.SetViewUp(start_vu)
1203    vp.camera.SetDistance(1065.168)
1204    vp.camera.SetClippingRange([517.354, 1757.322])
1205
1206    if include_faults:
1207        vp.show(
1208            msh,
1209            msh2,
1210            msh3,
1211            hist,
1212            brine_surface,
1213            oil_surface,
1214            gas_surface,
1215            fault_surface,
1216            text0,
1217            N=1,
1218            zoom=1.5,
1219        )
1220    else:
1221        vp.show(
1222            msh,
1223            msh2,
1224            msh3,
1225            hist,
1226            brine_surface,
1227            oil_surface,
1228            gas_surface,
1229            text0,
1230            N=1,
1231            zoom=1.5,
1232        )
1233
1234    # capture the scene from several vantage points
1235    delay = str(600 // camera_steps)
1236
1237    # define end camera position
1238    end_pos = [450.0, 275.0, 1150.0]
1239    end_vu = [-0.9, -0.3, 0.3]
1240
1241    n_images1 = 0
1242    camera_steps_1 = camera_steps // 3
1243    for icam in range(camera_steps):
1244        _scalar = float(icam) / (camera_steps - 1)
1245        _pos = np.array(start_pos) + (np.array(end_pos) - np.array(start_pos)) * _scalar
1246        _vu = np.array(start_vu) + (np.array(end_vu) - np.array(start_vu)) * _scalar
1247        vp.camera.SetPosition(_pos)
1248        vp.camera.SetViewUp(_vu)
1249        screenshot(
1250            filename=f"{cfg.temp_folder}/screenshot_{cfg.date_stamp}_{n_images1:03d}.png",
1251            scale=None,
1252            returnNumpy=False,
1253        )
1254        if cfg.verbose:
1255            print("\n" + str(n_images1), str(_pos), str(_vu))
1256        n_images1 += 1
1257
1258    start_pos = [450.0, 275.0, 1150.0]
1259    start_vu = [-0.9, -0.3, 0.3]
1260    end_pos = [275.0, 450.0, 1150.0]
1261    end_vu = [-0.3, -0.9, 0.3]
1262
1263    n_images2 = n_images1
1264    camera_steps_2 = camera_steps // 3
1265    for icam in range(camera_steps // 3):
1266        _scalar = float(icam) / (camera_steps - 1)
1267        _pos = np.array(start_pos) + (np.array(end_pos) - np.array(start_pos)) * _scalar
1268        _vu = np.array(start_vu) + (np.array(end_vu) - np.array(start_vu)) * _scalar
1269        vp.camera.SetPosition(_pos)
1270        vp.camera.SetViewUp(_vu)
1271        screenshot(
1272            filename=f"{cfg.temp_folder}/screenshot_{cfg.date_stamp}_{n_images2:03d}.png",
1273            scale=None,
1274            returnNumpy=False,
1275        )
1276        if cfg.verbose:
1277            print("\n" + str(n_images2), str(_pos), str(_vu))
1278        n_images2 += 1
1279
1280    start_pos = [275.0, 450.0, 1150.0]
1281    start_vu = [-0.3, -0.9, 0.3]
1282    end_vu = [-0.1, -0.4, 0.92]
1283    end_pos = [350.0, 825.0, 600.0]
1284
1285    n_images3 = n_images2
1286    camera_steps_3 = camera_steps // 3
1287    for icam in range(camera_steps // 3):
1288        _scalar = float(icam) / (camera_steps - 1)
1289        _pos = np.array(start_pos) + (np.array(end_pos) - np.array(start_pos)) * _scalar
1290        _vu = np.array(start_vu) + (np.array(end_vu) - np.array(start_vu)) * _scalar
1291        vp.camera.SetPosition(_pos)
1292        vp.camera.SetViewUp(_vu)
1293        screenshot(
1294            filename=f"{cfg.temp_folder}/screenshot_{cfg.date_stamp}_{n_images3:03d}.png",
1295            scale=None,
1296            returnNumpy=False,
1297        )
1298        if cfg.verbose:
1299            print("\n" + str(n_images3), str(_pos), str(_vu))
1300        n_images3 += 1
1301
1302    n_images4 = n_images3
1303    for icam in range(camera_steps // 6):
1304        fig_num = n_images3 + icam
1305        os.system(
1306            f"cp {cfg.temp_folder}/screenshot_{cfg.date_stamp}_{n_images3 - 1:03d}.png "
1307            + f"{cfg.temp_folder}/screenshot_{cfg.date_stamp}_{fig_num:03d}.png"
1308        )
1309        n_images4 += 1
1310
1311    # create animated gif from screenshots, convert requires ImageTools bin folder to be on $PATH
1312    output_imagename = os.path.join(
1313        cfg.work_subfolder, f"qc_image_3d_{cfg.date_stamp}.gif"
1314    )
1315    os.system(
1316        f"convert -delay {delay} {cfg.temp_folder}/screenshot_{cfg.date_stamp}_*.png -loop 0 \
1317          -dispose previous {output_imagename}"
1318    )
1319
1320    vp.close()
1321    vp.closeWindow()