datagenerator.Seismic

   1import os
   2import itertools
   3from multiprocessing import Pool
   4import numpy as np
   5
   6from tqdm import trange
   7from datagenerator.histogram_equalizer import normalize_seismic
   8from datagenerator.Geomodels import Geomodel
   9from datagenerator.util import write_data_to_hdf
  10from datagenerator.wavelets import generate_wavelet, plot_wavelets
  11from rockphysics.RockPropertyModels import select_rpm, RockProperties, EndMemberMixing
  12
  13
  14class SeismicVolume(Geomodel):
  15    """
  16    SeismicVolume Class
  17    -------------------
  18
  19    This class initializes the seismic volume that contains everything.
  20    """
  21    def __init__(self, parameters, faults, closures) -> None:
  22        """
  23        __init__
  24        --------
  25
  26        The initialization function.
  27
  28        Parameters
  29        ----------
  30        parameters : datagenerator.Parameters
  31            The parameters object of the project.
  32        faults : np.ndarray
  33            The faults array.
  34        closures : np.ndarray
  35            The closures array.
  36        
  37        Returns
  38        -------
  39        None
  40        """
  41        self.cfg = parameters
  42        self.faults = faults
  43        self.traps = closures
  44
  45        if self.cfg.model_qc_volumes:
  46            # Add 0 and 45 degree angles to the list of user-input angles
  47            self.angles = tuple(sorted(set((0,) + self.cfg.incident_angles + (45,))))
  48        else:
  49            self.angles = self.cfg.incident_angles
  50        if self.cfg.verbose:
  51            print(f"Angles: {self.angles}")
  52
  53        self.first_random_lyr = 20  # do not randomise shallow layers
  54
  55        self.rho = self.cfg.hdf_init(
  56            "rho", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  57        )
  58        self.vp = self.cfg.hdf_init(
  59            "vp", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  60        )
  61        self.vs = self.cfg.hdf_init(
  62            "vs", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  63        )
  64        self.rho_ff = self.cfg.hdf_init(
  65            "rho_ff", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  66        )
  67        self.vp_ff = self.cfg.hdf_init(
  68            "vp_ff", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  69        )
  70        self.vs_ff = self.cfg.hdf_init(
  71            "vs_ff", shape=self.cfg.h5file.root.ModelData.faulted_depth.shape
  72        )
  73
  74        seis_shape = (
  75            len(self.angles),
  76            *self.cfg.h5file.root.ModelData.faulted_depth.shape,
  77        )
  78        rfc_shape = (seis_shape[0], seis_shape[1], seis_shape[2], seis_shape[3] - 1)
  79        self.rfc_raw = self.cfg.hdf_init("rfc_raw", shape=rfc_shape)
  80        self.rfc_noise_added = self.cfg.hdf_init("rfc_noise_added", shape=rfc_shape)
  81
  82    def build_elastic_properties(self, mixing_method="inv_vel"):
  83        """
  84        Build Elastic Properties
  85        ------------------------
  86
  87        Creates Rho, Vp, Vs volumes using depth trends.
  88
  89        Parameters
  90        ----------
  91        mixing_method : str, optional
  92            Method for mixing Velocities. Defaults to "inv_vel".
  93
  94        """
  95        rpm = select_rpm(self.cfg)
  96        self.build_property_models_randomised_depth(
  97            rpm,
  98            mixing_method=mixing_method
  99        )
 100
 101    def bandlimit_volumes_wavelets(self, n_wavelets=10):
 102        # Pre-prepared n, m, f filter_specs from npy file
 103        fs_nr = np.load(self.cfg.wavelets[0])
 104        fs_md = np.load(self.cfg.wavelets[1])
 105        fs_fr = np.load(self.cfg.wavelets[2])
 106
 107        for wavelet_count in range(n_wavelets):
 108            low_freq_pctile = np.random.uniform(0, 100)
 109            mid_freq_pctile = np.random.uniform(0, 100)
 110            high_freq_pctile = np.random.uniform(0, 100)
 111            low_med_high_percentiles = (
 112                low_freq_pctile,
 113                mid_freq_pctile,
 114                high_freq_pctile,
 115            )
 116            print(f"low_med_high_percentiles: {low_med_high_percentiles}")
 117            f_nr, a_nr, w_nr = generate_wavelet(fs_nr, low_med_high_percentiles)
 118            f_md, a_md, w_md = generate_wavelet(fs_md, low_med_high_percentiles)
 119            f_fr, a_fr, w_fr = generate_wavelet(fs_fr, low_med_high_percentiles)
 120            _f = (f_nr, f_md, f_fr)
 121            _a = (a_nr, a_md, a_fr)
 122            _w = (w_nr, w_md, w_fr)
 123            labels = ["Near", "Mid", "Far"]
 124            pngname = os.path.join(
 125                self.cfg.work_subfolder, f"random_wavelets_{wavelet_count}.png"
 126            )
 127            plot_wavelets(_f, _a, _w, labels, pngname)
 128
 129            rfc_bandlimited_wavelet = np.zeros_like(self.rfc_noise_added)
 130            if self.cfg.model_qc_volumes:
 131                wavelets = [
 132                    w_nr,
 133                    w_nr,
 134                    w_md,
 135                    w_fr,
 136                    w_fr,
 137                ]  # use near for 0 and far for 45 deg cubes
 138            else:
 139                wavelets = [w_nr, w_md, w_fr]
 140
 141            for idx in range(rfc_bandlimited_wavelet.shape[0]):
 142                rfc_bandlimited_wavelet[idx, ...] = apply_wavelet(
 143                    self.rfc_noise_added[idx, ...], wavelets[idx]
 144                )
 145            if self.cfg.lateral_filter_size > 1:
 146                # Apply lateral smoothing filter
 147                rfc_bandlimited_wavelet = self.apply_lateral_filter(
 148                    rfc_bandlimited_wavelet
 149                )
 150            cumsum_wavelet_noise_added = self.apply_cumsum(rfc_bandlimited_wavelet)
 151            _ = self.write_final_cubes_to_disk(
 152                rfc_bandlimited_wavelet, f"seismicCubes_RFC_Wavelet_{wavelet_count}_"
 153            )
 154            _ = self.write_final_cubes_to_disk(
 155                cumsum_wavelet_noise_added,
 156                f"seismicCubes_cumsum_Wavelet_{wavelet_count}_",
 157            )
 158
 159    def postprocess_rfc_cubes(
 160        self, rfc_data, filename_suffix="", bb=False, stack=False
 161    ):
 162        """Postprocess the RFC cubes."""
 163        if bb:
 164            bandlimited = self.apply_bandlimits(rfc_data, (4.0, 90.0))
 165        else:
 166            bandlimited = self.apply_bandlimits(rfc_data)
 167        if self.cfg.lateral_filter_size > 1:
 168            bandlimited = self.apply_lateral_filter(bandlimited)
 169        if filename_suffix == "":
 170            fname = "seismicCubes_RFC_"
 171        else:
 172            fname = f"seismicCubes_RFC_{filename_suffix}_"
 173        _ = self.write_final_cubes_to_disk(bandlimited, fname)
 174        if stack:
 175            rfc_fullstack = self.stack_substacks(
 176                [bandlimited[x, ...] for x, _ in enumerate(self.angles)]
 177            )
 178            rfc_fullstack_scaled = self._scale_seismic(rfc_fullstack)
 179            self.write_cube_to_disk(
 180                rfc_fullstack_scaled,
 181                f"{fname}fullstack",
 182            )
 183
 184        cumsum = self.apply_cumsum(bandlimited)
 185        if filename_suffix == "":
 186            fname = "seismicCubes_cumsum_"
 187        else:
 188            fname = f"seismicCubes_cumsum_{filename_suffix}_"
 189        normalised_cumsum = self.write_final_cubes_to_disk(cumsum, fname)
 190        if stack:
 191            rai_fullstack = self.stack_substacks(
 192                [cumsum[x, ...] for x, _ in enumerate(self.angles)]
 193            )
 194            rai_fullstack_scaled = self._scale_seismic(rai_fullstack)
 195            self.write_cube_to_disk(
 196                rai_fullstack_scaled,
 197                f"{fname}fullstack",
 198            )
 199        return normalised_cumsum
 200
 201    def apply_augmentations(self, data, name="cumsum"):
 202        seabed = self.faults.faulted_depth_maps[..., 0] / self.cfg.digi
 203        augmented_data, augmented_labels = self.augment_data_and_labels(data, seabed)
 204        for i, ang in enumerate(self.cfg.incident_angles):
 205            data = augmented_data[i, ...]
 206            fname = f"seismicCubes_{name}_{ang}_degrees_normalized_augmented"
 207            self.write_cube_to_disk(data, fname)
 208            self.write_cube_to_disk(augmented_labels, "hc_closures_augmented")
 209
 210    def apply_rmo(self, data, name="cumsum_RMO"):
 211        """Apply RMO to the cumsum with noise. Note rmo application expects angle in last dimension"""
 212        rmo_indices = self.compute_randomRMO_1D(
 213            data.shape[-1], self.cfg.incident_angles, velocity_fraction=0.015
 214        )
 215        rmo_applied_data = np.zeros_like(np.moveaxis(data, 0, -1))
 216        for iline in range(self.cfg.cube_shape[0]):
 217            for xline in range(self.cfg.cube_shape[1]):
 218                rmo_applied_data[iline, xline, ...] = self.apply_rmo_augmentation(
 219                    np.moveaxis(data, 0, -1)[iline, xline, ...],
 220                    rmo_indices,
 221                    remove_mean_rmo=True,
 222                )
 223        rmo_applied_data = np.moveaxis(
 224            rmo_applied_data, -1, 0
 225        )  # Put the angles back into 1st dimension
 226        normalised_rmo_applied_data = self.write_final_cubes_to_disk(
 227            rmo_applied_data, f"seismicCubes_{name}_"
 228        )
 229        return normalised_rmo_applied_data
 230
 231    def build_seismic_volumes(self):
 232        """Build Seismic volumes.
 233
 234        * Use Zoeppritz to create RFC volumes for each required angle, add 0 & 45 degree stacks for QC
 235        * Add exponential weighted noise to the angle stacks
 236        * Apply bandpass filter
 237        * Apply lateral filter
 238        * Calculate cumsum
 239        * Create full stacks by stacking the user-input angles (i.e. not including 0, 45 degree stacks)
 240        * Write seismic volumes to disk
 241        """
 242        if self.cfg.verbose:
 243            print(
 244                f"Generating 3D Rock Properties Rho, Vp and Vs using depth trends for project {self.cfg.project}"
 245            )
 246        # Calculate Raw RFC from rho, vp, vs using Zoeppritz & write to HdF
 247        self.create_rfc_volumes()
 248        self.add_weighted_noise(self.faults.faulted_depth_maps)
 249        if hasattr(self.cfg, "wavelets"):
 250            self.bandlimit_volumes_wavelets(n_wavelets=1)
 251        normalised_cumsum = self.postprocess_rfc_cubes(
 252            self.rfc_noise_added[:], "", stack=True
 253        )
 254        self.apply_augmentations(normalised_cumsum, name="cumsum")
 255        if self.cfg.model_qc_volumes:
 256            _ = self.postprocess_rfc_cubes(self.rfc_raw[:], "noise_free", bb=False)
 257            if self.cfg.broadband_qc_volume:
 258                _ = self.postprocess_rfc_cubes(self.rfc_raw[:], "noise_free_bb", bb=True)
 259            normalised_cumsum_rmo = self.apply_rmo(normalised_cumsum)
 260            self.apply_augmentations(normalised_cumsum_rmo, name="cumsum_RMO")
 261
 262    def _scale_seismic(self, data):
 263        """Apply scaling factor to final seismic data volumes"""
 264        if data.ndim == 4:
 265            avg_stdev = np.mean([data[x, ...].std() for x in range(data.shape[0])])
 266        elif data.ndim == 3:
 267            avg_stdev = data.std()
 268        else:
 269            print(
 270                f"Input Data has {data.ndim} dimensions, should have 3 or 4 dimensions"
 271            )
 272        scaled_data = data * 100 / avg_stdev
 273        return scaled_data
 274
 275    def write_final_cubes_to_disk(self, dat, name):
 276        scaled_data = self._scale_seismic(dat)
 277
 278        # Apply scaling factors to N, M, F cubes from csv file before determining clip limits
 279        factor_near = self.cfg.rpm_scaling_factors["nearfactor"]
 280        factor_mid = self.cfg.rpm_scaling_factors["midfactor"]
 281        factor_far = self.cfg.rpm_scaling_factors["farfactor"]
 282        cube_incr = 0
 283        if self.cfg.model_qc_volumes and dat.shape[0] > 3:
 284            # 0 degrees has been added at start, and 45 at end
 285            cube_incr = 1
 286        scaled_data[0 + cube_incr, ...] *= factor_near
 287        scaled_data[1 + cube_incr, ...] *= factor_mid
 288        scaled_data[2 + cube_incr, ...] *= factor_far
 289
 290        for i, ang in enumerate(self.cfg.incident_angles):
 291            data = scaled_data[i, ...]
 292            fname = f"{name}_{ang}_degrees"
 293            self.write_cube_to_disk(data, fname)
 294        normed_data = normalize_seismic(scaled_data)
 295        for i, ang in enumerate(self.cfg.incident_angles):
 296            data = normed_data[i, ...]
 297            fname = f"{name}_{ang}_degrees_normalized"
 298            self.write_cube_to_disk(data, fname)
 299        return normed_data
 300
 301    @staticmethod
 302    def random_z_rho_vp_vs(dmin=-7, dmax=7):
 303        delta_z_rho = int(np.random.uniform(dmin, dmax))
 304        delta_z_vp = int(np.random.uniform(dmin, dmax))
 305        delta_z_vs = int(np.random.uniform(dmin, dmax))
 306        return delta_z_rho, delta_z_vp, delta_z_vs
 307
 308    def create_rfc_volumes(self):
 309        """
 310        Build a 3D array of PP-Reflectivity using Zoeppritz for each incident angle given.
 311
 312        :return: 4D array of PP-Reflectivity. 4th dimension holds reflectivities of each incident angle provided
 313        """
 314        theta = np.asanyarray(self.angles).reshape(
 315            (-1, 1)
 316        )  # Convert angles into a column vector
 317
 318        # Make empty array to store RPP via zoeppritz
 319        zoep = np.zeros(
 320            (self.rho.shape[0], self.rho.shape[1], self.rho.shape[2] - 1, theta.size),
 321            dtype="complex64",
 322        )
 323
 324        _rho = self.rho[:]
 325        _vp = self.vp[:]
 326        _vs = self.vs[:]
 327        if self.cfg.verbose:
 328            print("Checking for null values inside create_rfc_volumes:\n")
 329            print(f"Rho: {np.min(_rho)}, {np.max(_rho)}")
 330            print(f"Vp: {np.min(_vp)}, {np.max(_vp)}")
 331            print(f"Vs: {np.min(_vs)}, {np.max(_vs)}")
 332        # Doing this for each trace & all (5) angles is actually quicker than doing entire cubes for each angle
 333        for i in trange(
 334            self.rho.shape[0],
 335            desc=f"Calculating Zoeppritz for {len(self.angles)} angles",
 336        ):
 337            for j in range(self.rho.shape[1]):
 338                rho1, rho2 = _rho[i, j, :-1], _rho[i, j, 1:]
 339                vp1, vp2 = _vp[i, j, :-1], _vp[i, j, 1:]
 340                vs1, vs2 = _vs[i, j, :-1], _vs[i, j, 1:]
 341                rfc = RFC(vp1, vs1, rho1, vp2, vs2, rho2, theta)
 342                zoep[i, j, :, :] = rfc.zoeppritz_reflectivity().T
 343        # Set any voxels with imaginary parts to 0, since all transmitted energy travels along the reflector
 344        # zoep[np.where(np.imag(zoep) != 0)] = 0
 345        # discard complex values and set dtype to float64
 346        zoep = np.real(zoep).astype("float64")
 347        del _rho, _vp, _vs
 348
 349        # Move the angle from last dimension to first
 350        self.rfc_raw[:] = np.moveaxis(zoep, -1, 0)
 351
 352        if self.cfg.hdf_store:
 353            for n, d in zip(
 354                [
 355                    "qc_volume_rfc_raw_{}_degrees".format(
 356                        str(self.angles[x]).replace(".", "_")
 357                    )
 358                    for x in range(self.rfc_raw.shape[0])
 359                ],
 360                [self.rfc_raw[x, ...] for x in range(self.rfc_raw.shape[0])],
 361            ):
 362                write_data_to_hdf(n, d, self.cfg.hdf_master)
 363        if self.cfg.model_qc_volumes:
 364            # Write raw RFC values to disk
 365            _ = [
 366                self.write_cube_to_disk(
 367                    self.rfc_raw[x, ...],
 368                    f"qc_volume_rfc_raw_{self.angles[x]}_degrees",
 369                )
 370                for x in range(self.rfc_raw.shape[0])
 371            ]
 372            max_amp = np.max(np.abs(self.rfc_raw[:]))
 373            _ = [
 374                self.write_cube_to_disk(
 375                    self.rfc_raw[x, ...],
 376                    f"qc_volume_rfc_raw_{self.angles[x]}_degrees",
 377                )
 378                for x in range(self.rfc_raw.shape[0])
 379            ]
 380
 381    def noise_3d(self, cube_shape, verbose=False):
 382        if verbose:
 383            print("   ... inside noise3D")
 384
 385        noise_seed = np.random.randint(1, high=2 ** 32 - 1)
 386        sign_seed = np.random.randint(1, high=2 ** 32 - 1)
 387
 388        np.random.seed(noise_seed)
 389        noise3d = np.random.exponential(
 390            1.0 / 100.0, size=cube_shape[0] * cube_shape[1] * cube_shape[2]
 391        )
 392        np.random.seed(sign_seed)
 393        sign = np.random.binomial(
 394            1, 0.5, size=cube_shape[0] * cube_shape[1] * cube_shape[2]
 395        )
 396
 397        sign[sign == 0] = -1
 398        noise3d *= sign
 399        noise3d = noise3d.reshape((cube_shape[0], cube_shape[1], cube_shape[2]))
 400        return noise3d
 401
 402    def add_weighted_noise(self, depth_maps):
 403        """
 404        Apply 3D weighted random noise to ndarray.
 405
 406        :return: Sum of rfc_volumes and 3D noise model, weighted by angle
 407        """
 408
 409        from math import sin, cos
 410        from datagenerator.util import mute_above_seafloor
 411        from datagenerator.util import infill_surface
 412
 413        # Calculate standard deviation ratio to use based on user-input sn_db parameter
 414        if self.cfg.verbose:
 415            print(f"\n...adding random noise to {self.rfc_raw.shape[0]} cubes...")
 416            print(f"S/N ratio = {self.cfg.sn_db:.4f}")
 417        std_ratio = np.sqrt(10 ** (self.cfg.sn_db / 10.0))
 418
 419        # Compute standard deviation in bottom half of cube, and use this to scale the noise
 420        # if wb has holes (i.e. Nan or 0 values), fill holes using replacement with interpolated (NN) values
 421        wb_map = depth_maps[..., 0]
 422        # Add test for all 0's in wb_map in case z_shift has been applied
 423        if (
 424            np.isnan(np.min(wb_map)) or 0 in wb_map and not np.all(wb_map == 0)
 425        ):  # If wb contains nan's or at least one 0, this will return True
 426            wb_map = infill_surface(wb_map)
 427        wb_plus_15samples = wb_map / (self.cfg.digi + 15) * self.cfg.digi
 428
 429        # Use the Middle angle cube for normalisation
 430        if self.cfg.model_qc_volumes:
 431            # 0 and 45 degree cubes have been added
 432            norm_cube = self.rfc_raw[2, ...]
 433        else:
 434            norm_cube = self.rfc_raw[1, ...]
 435        data_below_seafloor = norm_cube[
 436            mute_above_seafloor(wb_plus_15samples, np.ones(norm_cube.shape, "float"))
 437            != 0.0
 438        ]
 439        data_std = data_below_seafloor.std()
 440
 441        # Create array to store the 3D noise model for each angle
 442        noise3d_cubes = np.zeros_like(self.rfc_raw[:])
 443
 444        # Add weighted stack of noise using Hilterman equation to each angle substack
 445        noise_0deg = self.noise_3d(self.rfc_raw.shape[1:], verbose=False)
 446        noise_45deg = self.noise_3d(self.rfc_raw.shape[1:], verbose=False)
 447        # Store the noise models
 448        self.noise_0deg = self.cfg.hdf_init("noise_0deg", shape=noise_0deg.shape)
 449        self.noise_45deg = self.cfg.hdf_init("noise_45deg", shape=noise_45deg.shape)
 450        self.noise_0deg[:] = noise_0deg
 451        self.noise_45deg[:] = noise_45deg
 452
 453        for x, ang in enumerate(self.angles):
 454            weighted_noise = noise_0deg * (cos(ang) ** 2) + noise_45deg * (
 455                sin(ang) ** 2
 456            )
 457            noise_std = weighted_noise.std()
 458            noise3d_cubes[x, ...] = weighted_noise * (data_std / noise_std) / std_ratio
 459            if self.cfg.verbose:
 460                print(
 461                    f"\t...Normalised noise3d for angle {ang}:\tMin: {noise3d_cubes[x, ...].min():.4f},"
 462                    f" mean: {noise3d_cubes[x, ...].mean():.4f}, max: {noise3d_cubes[x, ...].max():.4f},"
 463                    f" std: {noise3d_cubes[x, ...].std():.4f}"
 464                )
 465                print(
 466                    f"\tS/N ratio = {self.cfg.sn_db:3.1f} dB.\n\tstd_ratio = {std_ratio:.4f}"
 467                    f"\n\tdata_std = {data_std:.4f}\n\tnoise_std = {noise_std:.4f}"
 468                )
 469
 470        self.rfc_noise_added[:] = noise3d_cubes + self.rfc_raw[:]
 471
 472    def apply_bandlimits(self, data, frequencies=None):
 473        """
 474        Apply Butterworth Bandpass Filter to data
 475        :param data: 4D array of RFC values
 476        :param frequencies: Explicit frequency bounds (lower, upper)
 477        :return: 4D array of band-limited RFC
 478        """
 479        if self.cfg.verbose:
 480            print(f"Data Min: {np.min(data):.2f}, Data Max: {np.max(data):.2f}")
 481        dt = self.cfg.digi / 1000.0
 482        if (
 483            data.shape[-1] / self.cfg.infill_factor - self.cfg.pad_samples
 484            == self.cfg.cube_shape[-1]
 485        ):
 486            # Input data is infilled
 487            dt /= self.cfg.infill_factor
 488        if frequencies:  # if frequencies explicitly defined
 489            low = frequencies[0]
 490            high = frequencies[1]
 491        else:
 492            low = self.cfg.lowfreq
 493            high = self.cfg.highfreq
 494        b, a = derive_butterworth_bandpass(
 495            low, high, (dt * 1000.0), order=self.cfg.order
 496        )
 497        if self.cfg.verbose:
 498            print(f"\t... Low Frequency; {low:.2f} Hz, High Frequency: {high:.2f} Hz")
 499            print(
 500                f"\t... start_width: {1. / (dt * low):.4f}, end_width: {1. / (dt * high):.4f}"
 501            )
 502
 503        # Loop over 3D angle stack arrays
 504        if self.cfg.verbose:
 505            print(" ... shape of data being bandlimited = ", data.shape, "\n")
 506        num = data.shape[0]
 507        if self.cfg.verbose:
 508            print(f"Applying Bandpass to {num} cubes")
 509        if self.cfg.multiprocess_bp:
 510            # Much faster to use multiprocessing and apply band-limitation to entire cube at once
 511            with Pool(processes=min(num, os.cpu_count() - 2)) as p:
 512                iterator = zip(
 513                    [data[x, ...] for x in range(num)],
 514                    itertools.repeat(b),
 515                    itertools.repeat(a),
 516                )
 517                out_cubes_mp = p.starmap(self._run_bandpass_on_cubes, iterator)
 518            out_cubes = np.asarray(out_cubes_mp)
 519        else:
 520            # multiprocessing can fail using Python version 3.6 and very large arrays
 521            out_cubes = np.zeros_like(data)
 522            for idx in range(num):
 523                out_cubes[idx, ...] = self._run_bandpass_on_cubes(data[idx, ...], b, a)
 524
 525        return out_cubes
 526
 527    @staticmethod
 528    def _run_bandpass_on_cubes(data, b, a):
 529        out_cube = apply_butterworth_bandpass(data, b, a)
 530        return out_cube
 531
 532    def apply_lateral_filter(self, data):
 533        from scipy.ndimage import uniform_filter
 534
 535        n_filt = self.cfg.lateral_filter_size
 536        return uniform_filter(data, size=(0, n_filt, n_filt, 0))
 537
 538    def apply_cumsum(self, data):
 539        """Calculate cumulative sum on the Z-axis of input data
 540        Sipmap's cumsum also applies an Ormsby filter to the cumulatively summed data
 541        To replicate this, apply a butterworth filter to the data using 2, 100Hz frequency limits
 542        todo: make an ormsby function"""
 543        cumsum = data.cumsum(axis=-1)
 544        return self.apply_bandlimits(cumsum, (2.0, 100.0))
 545
 546    @staticmethod
 547    def stack_substacks(cube_list):
 548        """Stack cubes by taking mean"""
 549        return np.mean(cube_list, axis=0)
 550
 551    def augment_data_and_labels(self, normalised_seismic, seabed):
 552        from datagenerator.Augmentation import tz_stretch, uniform_stretch
 553
 554        hc_labels = self.cfg.h5file.root.ModelData.hc_labels[:]
 555        data, labels = tz_stretch(
 556            normalised_seismic,
 557            hc_labels[..., : self.cfg.cube_shape[2] + self.cfg.pad_samples - 1],
 558            seabed,
 559        )
 560        seismic, hc = uniform_stretch(data, labels)
 561        return seismic, hc
 562
 563    def compute_randomRMO_1D(
 564        self, nsamples, inc_angle_list, velocity_fraction=0.015, return_plot=False
 565    ):
 566        """
 567        compute 1D indices for to apply residual moveout via interpolation
 568        - nsamples is the array length in depth
 569        - inc_angle_list is a list of angle in degrees for angle stacks to which
 570        results will be applied
 571        - velocity_fraction is the maximum percent error in absolute value.
 572        for example, .03 would generate RMO based on velocity from .97 to 1.03 times correct velocity
 573        """
 574        from scipy.interpolate import UnivariateSpline
 575
 576        input_indices = np.arange(nsamples)
 577        if self.cfg.qc_plots:
 578            from datagenerator.util import import_matplotlib
 579
 580            plt = import_matplotlib()
 581            plt.close(1)
 582            plt.figure(1, figsize=(4.75, 6.0), dpi=150)
 583            plt.clf()
 584            plt.grid()
 585            plt.xlim((-10, 10))
 586            plt.ylim((nsamples, 0))
 587            plt.ylabel("depth (samples)")
 588            plt.xlabel("RMO (samples)")
 589            plt.savefig(
 590                os.path.join(self.cfg.work_subfolder, "rmo_1d.png"), format="png"
 591            )
 592        # generate 1 to 3 randomly located velocity fractions at random depth indices
 593        for _ in range(250):
 594            num_tie_points = np.random.randint(low=1, high=4)
 595            tie_fractions = np.random.uniform(
 596                low=-velocity_fraction, high=velocity_fraction, size=num_tie_points
 597            )
 598            if num_tie_points > 1:
 599                tie_indices = np.random.uniform(
 600                    low=0.25 * nsamples, high=0.75 * nsamples, size=num_tie_points
 601                ).astype("int")
 602
 603            else:
 604                tie_indices = int(
 605                    np.random.uniform(low=0.25 * nsamples, high=0.75 * nsamples)
 606                )
 607            tie_fractions = np.hstack(([0.0, 0.0], tie_fractions, [0.0, 0.0]))
 608            tie_indices = np.hstack(([0, 1], tie_indices, [nsamples - 2, nsamples - 1]))
 609            tie_indices = np.sort(tie_indices)
 610            if np.diff(tie_indices).min() <= 0.0:
 611                continue
 612            s = UnivariateSpline(tie_indices, tie_fractions, s=0.0)
 613            tie_fraction_array = s(input_indices)
 614            output_indices = np.zeros((nsamples, len(inc_angle_list)))
 615            for i, iangle in enumerate(inc_angle_list):
 616                output_indices[:, i] = input_indices * (
 617                    1.0
 618                    + tie_fraction_array * np.tan(float(iangle) * np.pi / 180.0) ** 2
 619                )
 620            if np.abs(output_indices[:, i] - input_indices).max() < 10.0:
 621                break
 622        for i, iangle in enumerate(inc_angle_list):
 623            if self.cfg.qc_plots:
 624                plt.plot(
 625                    output_indices[:, i] - input_indices,
 626                    input_indices,
 627                    label=str(inc_angle_list[i]),
 628                )
 629
 630                if i == len(inc_angle_list) - 1:
 631                    output_tie_indices = tie_indices * (
 632                        1.0
 633                        + tie_fraction_array[tie_indices]
 634                        * np.tan(float(iangle) * np.pi / 180.0) ** 2
 635                    )
 636                    plt.plot(output_tie_indices - tie_indices, tie_indices, "ko")
 637                    rmo_min = (output_indices[:, -1] - input_indices).min()
 638                    rmo_max = (output_indices[:, -1] - input_indices).max()
 639                    rmo_range = np.array((rmo_min, rmo_max)).round(1)
 640                    plt.title(
 641                        "simultated RMO arrival time offsets\n"
 642                        + "rmo range = "
 643                        + str(rmo_range),
 644                        fontsize=11,
 645                    )
 646                    plt.legend()
 647                    plt.tight_layout()
 648                plt.savefig(
 649                    os.path.join(self.cfg.work_subfolder, "rmo_arrival_time.png"),
 650                    format="png",
 651                )
 652
 653        if return_plot:
 654            return output_indices, plt
 655        else:
 656            return output_indices
 657
 658    def apply_rmo_augmentation(self, angle_gather, rmo_indices, remove_mean_rmo=False):
 659        """
 660        apply residual moveout via interpolation
 661        applied to 'angle_gather' with results from compute_randomRMO_1D in rmo_indices
 662        angle_gather and rmo_indices need to be 2D arrays with the same shape (n_samples x n_angles)
 663        - angle_gather is the array to which RMO is applied
 664        - rmo_indices is the array used to modify the two-way-times of angle_gather
 665        - remove_mean_rmo can be used to (when True) to apply relative RMO that keeps
 666        event energy at roughly the same TWT. for example to avoid applying RMO
 667        to both data and labels
 668        """
 669        if remove_mean_rmo:
 670            # apply de-trending to ensure that rmo does not change TWT
 671            # after rmo compared to before rmo 'augmentation'
 672            residual_times = rmo_indices.mean(axis=-1) - np.arange(
 673                angle_gather.shape[0]
 674            )
 675            residual_rmo_indices = rmo_indices * 1.0
 676            residual_rmo_indices[:, 0] -= residual_times
 677            residual_rmo_indices[:, 1] -= residual_times
 678            residual_rmo_indices[:, 2] -= residual_times
 679            _rmo_indices = residual_rmo_indices
 680        else:
 681            _rmo_indices = rmo_indices
 682
 683        rmo_angle_gather = np.zeros_like(angle_gather)
 684        for iangle in range(len(rmo_indices[-1])):
 685            rmo_angle_gather[:, iangle] = np.interp(
 686                range(angle_gather.shape[0]),
 687                _rmo_indices[:, iangle],
 688                angle_gather[:, iangle],
 689            )
 690
 691        return rmo_angle_gather
 692
 693    def build_property_models_randomised_depth(self, rpm, mixing_method="inv_vel"):
 694        """
 695        Build property models with randomised depth
 696        -------------------------------------------
 697        Builds property models with randomised depth.
 698
 699        v2 but with Backus moduli mixing instead of inverse velocity mixing
 700        mixing_method one of ["InverseVelocity", "BackusModuli"]
 701
 702        Parameters
 703        ----------
 704        rpm : str
 705            Rock properties model.
 706        mixing_method : str
 707            Mixing method for randomised depth.
 708        
 709        Returns
 710        -------
 711        None
 712        """
 713        layer_half_range = self.cfg.rpm_scaling_factors["layershiftsamples"]
 714        property_half_range = self.cfg.rpm_scaling_factors["RPshiftsamples"]
 715
 716        depth = self.cfg.h5file.root.ModelData.faulted_depth[:]
 717        lith = self.faults.faulted_lithology[:]
 718        net_to_gross = self.faults.faulted_net_to_gross[:]
 719
 720        oil_closures = self.cfg.h5file.root.ModelData.oil_closures[:]
 721        gas_closures = self.cfg.h5file.root.ModelData.gas_closures[:]
 722
 723        integer_faulted_age = (
 724            self.cfg.h5file.root.ModelData.faulted_age_volume[:] + 0.00001
 725        ).astype(int)
 726
 727        # Use empty class object to store all Rho, Vp, Vs volumes (randomised, fluid factor and non randomised)
 728        mixed_properties = ElasticProperties3D()
 729        mixed_properties.rho = np.zeros_like(lith)
 730        mixed_properties.vp = np.zeros_like(lith)
 731        mixed_properties.vs = np.zeros_like(lith)
 732        mixed_properties.sum_net = np.zeros_like(lith)
 733        cube_shape = lith.shape
 734
 735        if self.cfg.verbose:
 736            mixed_properties.rho_not_random = np.zeros_like(lith)
 737            mixed_properties.vp_not_random = np.zeros_like(lith)
 738            mixed_properties.vs_not_random = np.zeros_like(lith)
 739
 740        # water layer
 741        mixed_properties = self.water_properties(lith, mixed_properties)
 742
 743        mixed_properties.rho_ff = np.copy(mixed_properties.rho)
 744        mixed_properties.vp_ff = np.copy(mixed_properties.vp)
 745        mixed_properties.vs_ff = np.copy(mixed_properties.vs)
 746
 747        delta_z_lyrs = []
 748        delta_z_shales = []
 749        delta_z_brine_sands = []
 750        delta_z_oil_sands = []
 751        delta_z_gas_sands = []
 752
 753        for z in range(1, integer_faulted_age.max()):
 754            __i, __j, __k = np.where(integer_faulted_age == z)
 755
 756            if len(__k) > 0:
 757                if self.cfg.verbose:
 758                    print(
 759                        f"\n\n*********************\n ... layer number = {z}"
 760                        f"\n ... n2g (min,mean,max) = "
 761                        f"({np.min(net_to_gross[__i, __j, __k]).round(3)},"
 762                        f"{np.mean(net_to_gross[__i, __j, __k]).round(3)},"
 763                        f"{np.max(net_to_gross[__i, __j, __k]).round(3)},)"
 764                    )
 765
 766                delta_z_layer = self.get_delta_z_layer(z, layer_half_range, __k)
 767                delta_z_lyrs.append(delta_z_layer)
 768
 769            # shale required for all voxels
 770            i, j, k = np.where((net_to_gross >= 0.0) & (integer_faulted_age == z))
 771
 772            shale_voxel_count = len(k)
 773            delta_z_rho, delta_z_vp, delta_z_vs = self.get_delta_z_properties(
 774                z, property_half_range
 775            )
 776            delta_z_shales.append((delta_z_rho, delta_z_vp, delta_z_vs))
 777
 778            if len(k) > 0:
 779                _k_rho = list(
 780                    np.array(k + delta_z_layer + delta_z_rho).clip(
 781                        0, cube_shape[-1] - 10
 782                    )
 783                )
 784                _k_vp = list(
 785                    np.array(k + delta_z_layer + delta_z_vp).clip(
 786                        0, cube_shape[-1] - 10
 787                    )
 788                )
 789                _k_vs = list(
 790                    np.array(k + delta_z_layer + delta_z_vs).clip(
 791                        0, cube_shape[-1] - 10
 792                    )
 793                )
 794
 795                if self.cfg.verbose:
 796                    print(
 797                        f" ... shale: (delta_z_rho, delta_z_vp, delta_z_vs) = {delta_z_rho, delta_z_vp, delta_z_vs}"
 798                    )
 799                    print(f" ... shale: i = {np.mean(i)}")
 800                    print(f" ... shale: k = {np.mean(k)}")
 801                    print(f" ... shale: _k = {np.mean(_k_rho)}")
 802
 803                mixed_properties = self.calculate_shales(
 804                    xyz=(i, j, k),
 805                    shifts=(_k_rho, _k_vp, _k_vs),
 806                    props=mixed_properties,
 807                    rpm=rpm,
 808                    depth=depth,
 809                    z=z,
 810                )
 811
 812            # brine sand or mixture of brine sand and shale in same voxel
 813            i, j, k = np.where(
 814                (net_to_gross > 0.0)
 815                & (oil_closures == 0.0)
 816                & (gas_closures == 0.0)
 817                & (integer_faulted_age == z)
 818            )
 819            brine_voxel_count = len(k)
 820            delta_z_rho, delta_z_vp, delta_z_vs = self.get_delta_z_properties(
 821                z, property_half_range
 822            )
 823            delta_z_brine_sands.append((delta_z_rho, delta_z_vp, delta_z_vs))
 824
 825            if len(k) > 0:
 826                _k_rho = list(
 827                    np.array(k + delta_z_layer + delta_z_rho).clip(
 828                        0, cube_shape[-1] - 10
 829                    )
 830                )
 831                _k_vp = list(
 832                    np.array(k + delta_z_layer + delta_z_vp).clip(
 833                        0, cube_shape[-1] - 10
 834                    )
 835                )
 836                _k_vs = list(
 837                    np.array(k + delta_z_layer + delta_z_vs).clip(
 838                        0, cube_shape[-1] - 10
 839                    )
 840                )
 841
 842                if self.cfg.verbose:
 843                    print(
 844                        f"\n ... brine: (delta_z_rho, delta_z_vp, delta_z_vs) = {delta_z_rho, delta_z_vp, delta_z_vs}"
 845                    )
 846                    print(f" ... brine: i = {np.mean(i)}")
 847                    print(f" ... brine: k = {np.mean(k)}")
 848                    print(f" ... brine: _k = {np.mean(_k_rho)}")
 849
 850                mixed_properties = self.calculate_sands(
 851                    xyz=(i, j, k),
 852                    shifts=(_k_rho, _k_vp, _k_vs),
 853                    props=mixed_properties,
 854                    rpm=rpm,
 855                    depth=depth,
 856                    ng=net_to_gross,
 857                    z=z,
 858                    fluid="brine",
 859                    mix=mixing_method,
 860                )
 861
 862            # oil sands
 863            i, j, k = np.where(
 864                (net_to_gross > 0.0)
 865                & (oil_closures == 1.0)
 866                & (gas_closures == 0.0)
 867                & (integer_faulted_age == z)
 868            )
 869            oil_voxel_count = len(k)
 870            delta_z_rho, delta_z_vp, delta_z_vs = self.get_delta_z_properties(
 871                z, property_half_range
 872            )
 873            delta_z_oil_sands.append((delta_z_rho, delta_z_vp, delta_z_vs))
 874
 875            if len(k) > 0:
 876                _k_rho = list(
 877                    np.array(k + delta_z_layer + delta_z_rho).clip(
 878                        0, cube_shape[-1] - 10
 879                    )
 880                )
 881                _k_vp = list(
 882                    np.array(k + delta_z_layer + delta_z_vp).clip(
 883                        0, cube_shape[-1] - 10
 884                    )
 885                )
 886                _k_vs = list(
 887                    np.array(k + delta_z_layer + delta_z_vs).clip(
 888                        0, cube_shape[-1] - 10
 889                    )
 890                )
 891
 892                if self.cfg.verbose:
 893                    print(
 894                        "\n\n ... Perform fluid substitution for oil sands. ",
 895                        "\n ... Number oil voxels in closure = ",
 896                        mixed_properties.rho[i, j, k].shape,
 897                    )
 898                    print(
 899                        " ... closures_oil min/mean/max =  ",
 900                        oil_closures.min(),
 901                        oil_closures.mean(),
 902                        oil_closures.max(),
 903                    )
 904                    print(
 905                        "\n\n ... np.all(oil_sand_rho[:,:,:] == brine_sand_rho[:,:,:]) ",
 906                        np.all(
 907                            rpm.calc_oil_sand_properties(
 908                                depth[:], depth[:], depth[:]
 909                            ).rho
 910                            == rpm.calc_brine_sand_properties(
 911                                depth[:], depth[:], depth[:]
 912                            ).rho
 913                        ),
 914                    )
 915                    print(
 916                        " ... oil: (delta_z_rho, delta_z_vp, delta_z_vs) = "
 917                        + str((delta_z_rho, delta_z_vp, delta_z_vs))
 918                    )
 919
 920                if self.cfg.verbose:
 921                    print(" ... oil: i = " + str(np.mean(i)))
 922                    print(" ... oil: k = " + str(np.mean(k)))
 923                    print(" ... oil: _k = " + str(np.mean(_k_rho)))
 924
 925                mixed_properties = self.calculate_sands(
 926                    xyz=(i, j, k),
 927                    shifts=(_k_rho, _k_vp, _k_vs),
 928                    props=mixed_properties,
 929                    rpm=rpm,
 930                    depth=depth,
 931                    ng=net_to_gross,
 932                    z=z,
 933                    fluid="oil",
 934                    mix=mixing_method,
 935                )
 936
 937            # gas sands
 938            i, j, k = np.where(
 939                (net_to_gross > 0.0)
 940                & (oil_closures == 0.0)
 941                & (gas_closures == 1.0)
 942                & (integer_faulted_age == z)
 943            )
 944            gas_voxel_count = len(k)
 945            delta_z_rho, delta_z_vp, delta_z_vs = self.get_delta_z_properties(
 946                z, property_half_range
 947            )
 948            delta_z_gas_sands.append((delta_z_rho, delta_z_vp, delta_z_vs))
 949
 950            if len(k) > 0:
 951                _k_rho = list(
 952                    np.array(k + delta_z_layer + delta_z_rho).clip(
 953                        0, cube_shape[-1] - 10
 954                    )
 955                )
 956                _k_vp = list(
 957                    np.array(k + delta_z_layer + delta_z_vp).clip(
 958                        0, cube_shape[-1] - 10
 959                    )
 960                )
 961                _k_vs = list(
 962                    np.array(k + delta_z_layer + delta_z_vs).clip(
 963                        0, cube_shape[-1] - 10
 964                    )
 965                )
 966
 967                if self.cfg.verbose:
 968                    print(
 969                        "\n ... Perform fluid substitution for gas sands. ",
 970                        "\n ... Number gas voxels in closure = ",
 971                        mixed_properties.rho[i, j, k].shape,
 972                    )
 973                    print(
 974                        " ... closures_gas min/mean/max =  ",
 975                        gas_closures.min(),
 976                        gas_closures.mean(),
 977                        gas_closures.max(),
 978                    )
 979                    print(
 980                        " ... gas: (delta_z_rho, delta_z_vp, delta_z_vs) = "
 981                        + str((delta_z_rho, delta_z_vp, delta_z_vs))
 982                    )
 983
 984                if self.cfg.verbose:
 985                    print(" ... gas: i = " + str(np.mean(i)))
 986                    print(" ... gas: k = " + str(np.mean(k)))
 987                    print(" ... gas: _k = " + str(np.mean(_k_rho)))
 988
 989                mixed_properties = self.calculate_sands(
 990                    xyz=(i, j, k),
 991                    shifts=(_k_rho, _k_vp, _k_vs),
 992                    props=mixed_properties,
 993                    rpm=rpm,
 994                    depth=depth,
 995                    ng=net_to_gross,
 996                    z=z,
 997                    fluid="gas",
 998                    mix=mixing_method,
 999                )
1000
1001            # layer check
1002            __i, __j, __k = np.where(integer_faulted_age == z)
1003
1004            if len(__k) > 0:
1005                if self.cfg.verbose:
1006                    print(
1007                        "\n ... layer number = "
1008                        + str(z)
1009                        + "\n ... sum_net (min,mean,max) = "
1010                        + str(
1011                            (
1012                                mixed_properties.sum_net[__i, __j, __k].min(),
1013                                mixed_properties.sum_net[__i, __j, __k].mean(),
1014                                mixed_properties.sum_net[__i, __j, __k].max(),
1015                            )
1016                        )
1017                    )
1018                    print(
1019                        "\n ... layer number = "
1020                        + str(z)
1021                        + "\n ... (shale, brine, oil, gas) voxel_counts = "
1022                        + str(
1023                            (
1024                                shale_voxel_count,
1025                                brine_voxel_count,
1026                                oil_voxel_count,
1027                                gas_voxel_count,
1028                            )
1029                        )
1030                        + "\n ... shale_count = brine+oil+gas? "
1031                        + str(
1032                            shale_voxel_count
1033                            == brine_voxel_count + oil_voxel_count + gas_voxel_count
1034                        )
1035                        + "\n*********************"
1036                    )
1037
1038        # overwrite rho, vp, vs for salt if required
1039        if self.cfg.include_salt:
1040            # Salt. Set lith = 2.0
1041            i, j, k = np.where(lith == 2.0)
1042            mixed_properties.rho[i, j, k] = 2.17  # divide by 2 since lith = 2.0
1043            mixed_properties.vp[i, j, k] = 4500.0
1044            mixed_properties.vs[i, j, k] = 2250.0
1045            # Include fluidfactor rho, vp, vs inside salt body
1046            mixed_properties.rho_ff[i, j, k] = mixed_properties.rho[i, j, k]
1047            mixed_properties.vp_ff[i, j, k] = mixed_properties.vp[i, j, k]
1048            mixed_properties.vs_ff[i, j, k] = mixed_properties.vs[i, j, k]
1049
1050        if self.cfg.verbose:
1051            print("Checking for null values inside build_randomised_properties:\n")
1052            print(
1053                f"Rho: {np.min(mixed_properties.rho)}, {np.max(mixed_properties.rho)}"
1054            )
1055            print(f"Vp: {np.min(mixed_properties.vp)}, {np.max(mixed_properties.vp)}")
1056            print(f"Vs: {np.min(mixed_properties.vs)}, {np.max(mixed_properties.vs)}")
1057
1058        # Fix empty voxels at base of property volumes
1059        mixed_properties = self.fix_zero_values_at_base(mixed_properties)
1060
1061        if self.cfg.verbose:
1062            print(
1063                "Checking for null values inside build_randomised_properties AFTER fix:\n"
1064            )
1065            print(
1066                f"Rho: {np.min(mixed_properties.rho)}, {np.max(mixed_properties.rho)}"
1067            )
1068            print(f"Vp: {np.min(mixed_properties.vp)}, {np.max(mixed_properties.vp)}")
1069            print(f"Vs: {np.min(mixed_properties.vs)}, {np.max(mixed_properties.vs)}")
1070
1071        mixed_properties = self.apply_scaling_factors(
1072            mixed_properties, net_to_gross, lith
1073        )
1074        self.rho[:] = mixed_properties.rho
1075        self.vp[:] = mixed_properties.vp
1076        self.vs[:] = mixed_properties.vs
1077
1078        self.rho_ff[:] = mixed_properties.rho_ff
1079        self.vp_ff[:] = mixed_properties.vp_ff
1080        self.vs_ff[:] = mixed_properties.vs_ff
1081
1082        if self.cfg.qc_plots:
1083            print("\nCreating RPM qc plots")
1084            from rockphysics.RockPropertyModels import rpm_qc_plots
1085
1086            rpm_qc_plots(self.cfg, rpm)
1087
1088        if self.cfg.model_qc_volumes:
1089            self.write_property_volumes_to_disk()
1090
1091    def water_properties(self, lith, properties):
1092        i, j, k = np.where(lith < 0.0)
1093        properties.rho[i, j, k] = 1.028
1094        properties.vp[i, j, k] = 1500.0
1095        properties.vs[i, j, k] = 1000.0
1096        return properties
1097
1098    def get_delta_z_layer(self, z, half_range, z_cells):
1099        if z > self.first_random_lyr:
1100            delta_z_layer = int(np.random.uniform(-half_range, half_range))
1101        else:
1102            delta_z_layer = 0
1103        if self.cfg.verbose:
1104            print(f" .... Layer {z}: voxel_count = {len(z_cells)}")
1105            print(f" .... Layer {z}: delta_z_layer = {delta_z_layer}")
1106            print(
1107                f" .... Layer {z}: z-range (m): {np.min(z_cells) * self.cfg.digi}, "
1108                f"{np.max(z_cells) * self.cfg.digi}"
1109            )
1110        return delta_z_layer
1111
1112    def get_delta_z_properties(self, z, half_range):
1113        if z > self.first_random_lyr:
1114            delta_z_rho, delta_z_vp, delta_z_vs = self.random_z_rho_vp_vs(
1115                dmin=-half_range, dmax=half_range
1116            )
1117        else:
1118            delta_z_rho, delta_z_vp, delta_z_vs = (0, 0, 0)
1119        return delta_z_rho, delta_z_vp, delta_z_vs
1120
1121    def calculate_shales(self, xyz, shifts, props, rpm, depth, z):
1122        # Shales required for all voxels, other than water
1123        # Calculate the properties, select how to mix with other facies later
1124        i, j, k = xyz
1125        k_rho, k_vp, k_vs = shifts
1126
1127        z_rho = depth[i, j, k_rho]
1128        z_vp = depth[i, j, k_vp]
1129        z_vs = depth[i, j, k_vs]
1130        shales = rpm.calc_shale_properties(z_rho, z_vp, z_vs)
1131
1132        # Do not use net to gross here.
1133        # Since every voxel will have some shale in, apply the net to gross
1134        # when combining shales and sands
1135        # _ng = (1.0 - ng[i, j, k])
1136
1137        props.rho[i, j, k] = shales.rho  # * _ng
1138        props.rho_ff[i, j, k] = shales.rho  # * _ng
1139        props.vp[i, j, k] = shales.vp  # * _ng
1140        props.vp_ff[i, j, k] = shales.vp  # * _ng
1141        props.vs[i, j, k] = shales.vs  # * _ng
1142        props.vs_ff[i, j, k] = shales.vs  # * _ng
1143        # props.sum_net[i, j, k] = _ng
1144
1145        if self.cfg.verbose and z > self.first_random_lyr:
1146            # Calculate non-randomised properties and differences
1147            _z0 = depth[i, j, k]
1148            _shales0 = rpm.calc_shale_properties(_z0, _z0, _z0)
1149            props.rho_not_random[i, j, k] = _shales0.rho  # * _ng
1150            props.vp_not_random[i, j, k] = _shales0.vp  # * _ng
1151            props.vs_not_random[i, j, k] = _shales0.vs  # * _ng
1152
1153            delta_rho = 1.0 - (props.rho[i, j, k] / props.rho_not_random[i, j, k])
1154            delta_vp = 1.0 - (props.vp[i, j, k] / props.vp_not_random[i, j, k])
1155            delta_vs = 1.0 - (props.vs[i, j, k] / props.vs_not_random[i, j, k])
1156            pct_change_rho = delta_rho.mean().round(3)
1157            pct_change_vp = delta_vp.mean().round(3)
1158            pct_change_vs = delta_vs.mean().round(3)
1159            print(
1160                f" ... shale: randomization percent change (rho,vp,vs) = "
1161                f"{pct_change_rho}, {pct_change_vp }, {pct_change_vs}"
1162            )
1163            print(
1164                f" ... shale: min/max pre-randomization (rho, vp, vs)  = "
1165                f"{np.min(props.rho_not_random[i, j, k]):.3f} - "
1166                f"{np.max(props.rho_not_random[i, j, k]):.3f}, "
1167                f"{np.min(props.vp_not_random[i, j, k]):.2f} - "
1168                f"{np.max(props.vp_not_random[i, j, k]):.2f}, "
1169                f"{np.min(props.vs_not_random[i, j, k]):.2f} - "
1170                f"{np.max(props.vs_not_random[i, j, k]):.2f}"
1171            )
1172            print(
1173                f" ... shale: min/max post-randomization (rho, vp, vs) = "
1174                f"{np.min(props.rho[i, j, k]):.3f} - "
1175                f"{np.max(props.rho[i, j, k]):.3f}, "
1176                f"{np.min(props.vp[i, j, k]):.2f} - "
1177                f"{np.max(props.vp[i, j, k]):.2f}, "
1178                f"{np.min(props.vs[i, j, k]):.2f} - "
1179                f"{np.max(props.vs[i, j, k]):.2f}"
1180            )
1181
1182        return props
1183
1184    def calculate_sands(
1185        self, xyz, shifts, props, rpm, depth, ng, z, fluid, mix="inv_vel"
1186    ):
1187        # brine sand or mixture of brine sand and shale in same voxel
1188        # - perform velocity sums in slowness or via backus moduli mixing
1189        i, j, k = xyz
1190        k_rho, k_vp, k_vs = shifts
1191
1192        z_rho = depth[i, j, k_rho]
1193        z_vp = depth[i, j, k_vp]
1194        z_vs = depth[i, j, k_vs]
1195        if fluid == "brine":
1196            sands = rpm.calc_brine_sand_properties(z_rho, z_vp, z_vs)
1197        elif fluid == "oil":
1198            sands = rpm.calc_oil_sand_properties(z_rho, z_vp, z_vs)
1199        elif fluid == "gas":
1200            sands = rpm.calc_gas_sand_properties(z_rho, z_vp, z_vs)
1201
1202        shales = RockProperties(
1203            rho=props.rho[i, j, k], vp=props.vp[i, j, k], vs=props.vs[i, j, k]
1204        )
1205        rpmix = EndMemberMixing(shales, sands, ng[i, j, k])
1206
1207        if mix == "inv_vel":
1208            rpmix.inverse_velocity_mixing()
1209        else:
1210            rpmix.backus_moduli_mixing()
1211
1212        props.sum_net[i, j, k] += ng[i, j, k]
1213        props.rho[i, j, k] = rpmix.rho
1214        props.vp[i, j, k] = rpmix.vp
1215        props.vs[i, j, k] = rpmix.vs
1216
1217        if self.cfg.verbose and z > self.first_random_lyr:
1218            # Calculate non-randomised properties and differences
1219            _z0 = depth[i, j, k]
1220            if fluid == "brine":
1221                _sands0 = rpm.calc_brine_sand_properties(_z0, _z0, _z0)
1222            elif fluid == "oil":
1223                _sands0 = rpm.calc_oil_sand_properties(_z0, _z0, _z0)
1224            elif fluid == "gas":
1225                _sands0 = rpm.calc_gas_sand_properties(_z0, _z0, _z0)
1226
1227            rpmix_0 = EndMemberMixing(shales, _sands0, ng[i, j, k])
1228            if mix == "inv_vel":
1229                # Apply Inverse Velocity mixing
1230                rpmix_0.inverse_velocity_mixing()
1231            else:
1232                # Apply Backus Moduli mixing
1233                rpmix_0.backus_moduli_mixing()
1234            props.rho_not_random[i, j, k] = rpmix_0.rho
1235            props.vp_not_random[i, j, k] = rpmix_0.vp
1236            props.vs_not_random[i, j, k] = rpmix_0.vs
1237
1238            delta_rho = 1.0 - (props.rho[i, j, k] / props.rho_not_random[i, j, k])
1239            delta_vp = 1.0 - (props.vp[i, j, k] / props.vp_not_random[i, j, k])
1240            delta_vs = 1.0 - (props.vs[i, j, k] / props.vs_not_random[i, j, k])
1241            pct_change_rho = delta_rho.mean().round(3)
1242            pct_change_vp = delta_vp.mean().round(3)
1243            pct_change_vs = delta_vs.mean().round(3)
1244            print(
1245                f" ... {fluid}: randomization percent change (rho,vp,vs) = "
1246                f"{pct_change_rho}, {pct_change_vp }, {pct_change_vs}"
1247            )
1248            print(
1249                f" ... {fluid}: min/max pre-randomization (rho, vp, vs)  = "
1250                f"{np.min(props.rho_not_random[i, j, k]):.3f} - "
1251                f"{np.max(props.rho_not_random[i, j, k]):.3f}, "
1252                f"{np.min(props.vp_not_random[i, j, k]):.2f} - "
1253                f"{np.max(props.vp_not_random[i, j, k]):.2f}, "
1254                f"{np.min(props.vs_not_random[i, j, k]):.2f} - "
1255                f"{np.max(props.vs_not_random[i, j, k]):.2f}"
1256            )
1257            print(
1258                f" ... {fluid}: min/max post-randomization (rho, vp, vs) = "
1259                f"{np.min(props.rho[i, j, k]):.3f} - "
1260                f"{np.max(props.rho[i, j, k]):.3f}, "
1261                f"{np.min(props.vp[i, j, k]):.2f} - "
1262                f"{np.max(props.vp[i, j, k]):.2f}, "
1263                f"{np.min(props.vs[i, j, k]):.2f} - "
1264                f"{np.max(props.vs[i, j, k]):.2f}"
1265            )
1266
1267        return props
1268
1269    @staticmethod
1270    def fix_zero_values_at_base(props):
1271        """Check for zero values at base of property volumes and replace with
1272        shallower values if present
1273
1274        Args:
1275            a ([type]): [description]
1276            b ([type]): [description]
1277            c ([type]): [description]
1278        """
1279        for vol in [
1280            props.rho,
1281            props.vp,
1282            props.vs,
1283            props.rho_ff,
1284            props.vp_ff,
1285            props.vs_ff,
1286        ]:
1287            for (x, y, z) in np.argwhere(vol == 0.0):
1288                vol[x, y, z] = vol[x, y, z - 1]
1289        return props
1290
1291    def apply_scaling_factors(self, props, ng, lith):
1292        """Apply final random scaling factors to sands and shales"""
1293        rho_factor_shale = self.cfg.rpm_scaling_factors["shalerho_factor"]
1294        vp_factor_shale = self.cfg.rpm_scaling_factors["shalevp_factor"]
1295        vs_factor_shale = self.cfg.rpm_scaling_factors["shalevs_factor"]
1296        rho_factor_sand = self.cfg.rpm_scaling_factors["sandrho_factor"]
1297        vp_factor_sand = self.cfg.rpm_scaling_factors["sandvp_factor"]
1298        vs_factor_sand = self.cfg.rpm_scaling_factors["sandvs_factor"]
1299
1300        if self.cfg.verbose:
1301            print(
1302                f"\n... Additional random scaling factors: "
1303                f"\n ... Shale Rho {rho_factor_shale:.3f}, Vp {vp_factor_shale:.3f}, Vs {vs_factor_shale:.3f}"
1304                f"\n ... Sand Rho {rho_factor_sand:.3f}, Vp {vp_factor_sand:.3f}, Vs {vs_factor_sand:.3f}"
1305            )
1306
1307        # Apply final scaling factors to shale/sand properties
1308        props.rho[(ng <= 1.0e-2) & (lith < 2.0)] = (
1309            props.rho[(ng <= 1.0e-2) & (lith < 2.0)] * rho_factor_shale
1310        )
1311        props.vp[(ng <= 1.0e-2) & (lith < 2.0)] = (
1312            props.vp[(ng <= 1.0e-2) & (lith < 2.0)] * vp_factor_shale
1313        )
1314        props.vs[(ng <= 1.0e-2) & (lith < 2.0)] = (
1315            props.vs[(ng <= 1.0e-2) & (lith < 2.0)] * vs_factor_shale
1316        )
1317        props.rho[(ng > 1.0e-2) & (lith < 2.0)] = (
1318            props.rho[(ng > 1.0e-2) & (lith < 2.0)] * rho_factor_sand
1319        )
1320        props.vp[(ng > 1.0e-2) & (lith < 2.0)] = (
1321            props.vp[(ng > 1.0e-2) & (lith < 2.0)] * vp_factor_sand
1322        )
1323        props.vs[(ng > 1.0e-2) & (lith < 2.0)] = (
1324            props.vs[(ng > 1.0e-2) & (lith < 2.0)] * vs_factor_sand
1325        )
1326        # Apply same factors to the fluid-factor property cubes
1327        props.rho_ff[(ng <= 1.0e-2) & (lith < 2.0)] = (
1328            props.rho_ff[(ng <= 1.0e-2) & (lith < 2.0)] * rho_factor_shale
1329        )
1330        props.vp_ff[(ng <= 1.0e-2) & (lith < 2.0)] = (
1331            props.vp_ff[(ng <= 1.0e-2) & (lith < 2.0)] * vp_factor_shale
1332        )
1333        props.vs_ff[(ng <= 1.0e-2) & (lith < 2.0)] = (
1334            props.vs_ff[(ng <= 1.0e-2) & (lith < 2.0)] * vs_factor_shale
1335        )
1336        props.rho_ff[(ng > 1.0e-2) & (lith < 2.0)] = (
1337            props.rho_ff[(ng > 1.0e-2) & (lith < 2.0)] * rho_factor_sand
1338        )
1339        props.vp_ff[(ng > 1.0e-2) & (lith < 2.0)] = (
1340            props.vp_ff[(ng > 1.0e-2) & (lith < 2.0)] * vp_factor_sand
1341        )
1342        props.vs_ff[(ng > 1.0e-2) & (lith < 2.0)] = (
1343            props.vs_ff[(ng > 1.0e-2) & (lith < 2.0)] * vs_factor_sand
1344        )
1345        return props
1346
1347    def write_property_volumes_to_disk(self):
1348        """Write Rho, Vp, Vs volumes to disk."""
1349        self.write_cube_to_disk(
1350            self.rho[:],
1351            "qc_volume_rho",
1352        )
1353        self.write_cube_to_disk(
1354            self.vp[:],
1355            "qc_volume_vp",
1356        )
1357        self.write_cube_to_disk(
1358            self.vs[:],
1359            "qc_volume_vs",
1360        )
1361
1362
1363class ElasticProperties3D:
1364    """Empty class to hold 3D property cubes Rho, Vp, Vs, Rho_ff, Vp_ff, Vs_ff and sum_net"""
1365
1366    def __init__(self):
1367        self.rho = None
1368        self.vp = None
1369        self.vs = None
1370        self.rho_ff = None
1371        self.vp_ff = None
1372        self.vs_ff = None
1373        self.sum_net = None
1374        self.rho_not_random = None
1375        self.vp_not_random = None
1376        self.vs_not_random = None
1377
1378
1379class RFC:
1380    """
1381    Reflection Coefficient object
1382    Contains methods for calculating the reflection coefficient at an interface
1383    from input Vp, Vs, Rho and incident angle
1384    Angle should be given in degrees and is converted to radians internally
1385    """
1386
1387    def __init__(
1388        self, vp_upper, vs_upper, rho_upper, vp_lower, vs_lower, rho_lower, theta
1389    ):
1390        self.theta = theta
1391        self.vp_upper = vp_upper
1392        self.vs_upper = vs_upper
1393        self.rho_upper = rho_upper
1394        self.vp_lower = vp_lower
1395        self.vs_lower = vs_lower
1396        self.rho_lower = rho_lower
1397
1398    def normal_incidence_rfc(self):
1399        z0 = self.rho_upper * self.vp_upper
1400        z1 = self.rho_lower * self.vp_lower
1401        return (z1 - z0) / (z1 + z0)
1402
1403    def shuey(self):
1404        """Use 3-term Shuey to calculate RFC."""
1405        radians = self._degrees_to_radians()  # Angle to radians
1406        sin_squared = np.sin(radians) ** 2
1407        tan_squared = np.tan(radians) ** 2
1408
1409        # Delta Vp, Vs, Rho
1410        d_vp = self.vp_lower - self.vp_upper
1411        d_vs = self.vs_lower - self.vs_upper
1412        d_rho = self.rho_lower - self.rho_upper
1413        avg_vp = np.mean([self.vp_lower, self.vp_upper])
1414        avg_vs = np.mean([self.vs_lower, self.vs_upper])
1415        avg_rho = np.mean([self.rho_lower, self.rho_upper])
1416
1417        # 3 Term Shuey
1418        r0 = 0.5 * (d_vp / avg_vp + d_rho / avg_rho)
1419        g = 0.5 * (d_vp / avg_vp) - 2.0 * avg_vs ** 2 / avg_vp ** 2 * (
1420            d_rho / avg_rho + 2.0 * d_vs / avg_vs
1421        )
1422        f = 0.5 * (d_vp / avg_vp)
1423
1424        return r0 + g * sin_squared + f * (tan_squared - sin_squared)
1425
1426    def zoeppritz_reflectivity(self):
1427        """Calculate PP reflectivity using Zoeppritz equation"""
1428        theta = self._degrees_to_radians(
1429            dtype="complex"
1430        )  # Convert angle to radians, as complex value
1431        p = np.sin(theta) / self.vp_upper  # ray param
1432        theta2 = np.arcsin(p * self.vp_lower)  # Transmission angle of P-wave
1433        phi1 = np.arcsin(p * self.vs_upper)  # Reflection angle of converted S-wave
1434        phi2 = np.arcsin(p * self.vs_lower)  # Transmission angle of converted S-wave
1435
1436        a = self.rho_lower * (1 - 2 * np.sin(phi2) ** 2.0) - self.rho_upper * (
1437            1 - 2 * np.sin(phi1) ** 2.0
1438        )
1439        b = (
1440            self.rho_lower * (1 - 2 * np.sin(phi2) ** 2.0)
1441            + 2 * self.rho_upper * np.sin(phi1) ** 2.0
1442        )
1443        c = (
1444            self.rho_upper * (1 - 2 * np.sin(phi1) ** 2.0)
1445            + 2 * self.rho_lower * np.sin(phi2) ** 2.0
1446        )
1447        d = 2 * (
1448            self.rho_lower * self.vs_lower ** 2 - self.rho_upper * self.vs_upper ** 2
1449        )
1450
1451        e = (b * np.cos(theta) / self.vp_upper) + (c * np.cos(theta2) / self.vp_lower)
1452        f = (b * np.cos(phi1) / self.vs_upper) + (c * np.cos(phi2) / self.vs_lower)
1453        g = a - d * np.cos(theta) / self.vp_upper * np.cos(phi2) / self.vs_lower
1454        h = a - d * np.cos(theta2) / self.vp_lower * np.cos(phi1) / self.vs_upper
1455
1456        d = e * f + g * h * p ** 2
1457
1458        zoep_pp = (
1459            f
1460            * (
1461                b * (np.cos(theta) / self.vp_upper)
1462                - c * (np.cos(theta2) / self.vp_lower)
1463            )
1464            - h
1465            * p ** 2
1466            * (a + d * (np.cos(theta) / self.vp_upper) * (np.cos(phi2) / self.vs_lower))
1467        ) * (1 / d)
1468
1469        return np.squeeze(zoep_pp)
1470
1471    def _degrees_to_radians(self, dtype="float"):
1472        """Convert angle in degrees to ange in radians. If dtype != 'float', return a complex dtype array"""
1473        if dtype == "float":
1474            return np.radians(self.theta)
1475        else:
1476            return np.radians(self.theta).astype(np.complex)
1477
1478
1479def derive_butterworth_bandpass(lowcut, highcut, digitisation, order=4):
1480    from scipy.signal import butter
1481
1482    fs = 1.0 / (digitisation / 1000.0)
1483    nyq = 0.5 * fs
1484    low = lowcut / nyq
1485    high = highcut / nyq
1486    b, a = butter(order, [low, high], btype="bandpass", output="ba")
1487    return b, a
1488
1489
1490def apply_butterworth_bandpass(data, b, a):
1491    """Apply Butterworth bandpass to data"""
1492    from scipy.signal import tf2zpk, filtfilt
1493
1494    # Use irlen to remove artefacts generated at base of cubes during bandlimitation
1495    _, p, _ = tf2zpk(b, a)
1496    eps = 1e-9
1497    r = np.max(np.abs(p))
1498    approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r)))
1499    y = filtfilt(b, a, data, method="gust", irlen=approx_impulse_len)
1500    # y = filtfilt(b, a, data)
1501    return y
1502
1503
1504def trim_triangular(low, mid, high):
1505    ###
1506    ### use trimmed triangular distributions
1507    ###
1508    from numpy.random import triangular
1509
1510    for _ in range(50):
1511        num = triangular(2 * low - mid, mid, 2 * high - mid)
1512        if low <= num <= high:
1513            break
1514    return num
1515
1516
1517def apply_wavelet(cube, wavelet):
1518    filtered_cube = np.zeros_like(cube)
1519    for i in range(cube.shape[0]):
1520        for j in range(cube.shape[1]):
1521            filtered_cube[i, j, :] = np.convolve(cube[i, j, :], wavelet, mode="same")
1522    return filtered_cube