datagenerator.Geomodels 

  1import os
  2import numpy as np
  3from datagenerator.Parameters import Parameters
  4from datagenerator.util import next_odd
  5from scipy.ndimage import maximum_filter
  6
  7
  8class Geomodel:
  9    """
 10        Geomodel
 11        --------------------------
 12        The class of the Geomodel object.
 13
 14        This class contains all the items that make up the Geologic model.
 15
 16        Parameters
 17        ----------
 18        parameters : datagenerator.Parameters
 19            Parameter object storing all model parameters.
 20        depth_maps : np.ndarray
 21            A numpy array containing the depth maps.
 22        onlap_horizon_list : list
 23            A list of the onlap horizons.
 24        facies : np.ndarray
 25            The generated facies.
 26
 27        Returns
 28        -------
 29        None
 30    """
 31    def __init__(
 32        self,
 33        parameters: Parameters,
 34        depth_maps: np.ndarray,
 35        onlap_horizon_list: list,
 36        facies: np.ndarray
 37    ) -> None:
 38        """__init__
 39
 40        Initializer for the Geomodel class.
 41
 42        Parameters
 43        ----------
 44        parameters : datagenerator.Parameters
 45            Parameter object storing all model parameters.
 46        depth_maps : np.ndarray
 47            A numpy array containing the depth maps.
 48        onlap_horizon_list : list
 49            A list of the onlap horizons.
 50        facies : np.ndarray
 51            The generated facies.
 52        """
 53        self.cfg = parameters
 54        self.depth_maps = depth_maps
 55        self.onlap_horizon_list = onlap_horizon_list
 56        self.facies = facies
 57
 58        # Initialise geomodels
 59        cube_shape = (
 60            self.cfg.cube_shape[0],
 61            self.cfg.cube_shape[1],
 62            self.cfg.cube_shape[2] + self.cfg.pad_samples,
 63        )
 64        self.geologic_age = self.cfg.hdf_init("geologic_age_prefault", shape=cube_shape)
 65        self.onlap_segments = self.cfg.hdf_init(
 66            "onlap_segments_prefault", shape=cube_shape
 67        )
 68        self.faulted_lithology = self.cfg.hdf_init(
 69            "lithology_prefault", shape=cube_shape
 70        )
 71        self.geomodel_ng = self.cfg.hdf_init("net_to_gross_prefault", shape=cube_shape)
 72        self.faulted_depth = self.cfg.hdf_init("depth_prefault", shape=cube_shape)
 73        self.faulted_depth_randomised = self.cfg.hdf_init(
 74            "depth_randomised_prefault", shape=cube_shape
 75        )
 76
 77        # Channel volumes
 78        if self.cfg.include_channels:
 79            self.floodplain_shale = None
 80            self.channel_fill = None
 81            self.shale_channel_drape = None
 82            self.levee = None
 83            self.crevasse = None
 84            self.channel_segments = None
 85
 86    def build_unfaulted_geomodels(self):
 87        """
 88            Build unfaulted geomodels.
 89            --------------------------
 90            A method that builds unfaulted geomodels.
 91
 92            This method does the following:
 93
 94            * Build geologic age cube from depth horizons
 95            * Build onlap segmentation cube
 96            * If channels are turned on, use fluvsim fortran code to build:
 97            - floodplain shale cube
 98            - channel fill cube
 99            - shale channel drape
100            - levee cube
101            - crevasse cube
102
103            Parameters
104            ----------
105            None
106
107            Returns
108            -------
109            None
110        """
111        self.geologic_age[:] = self.create_geologic_age_3d_from_infilled_horizons(
112            self.depth_maps
113        )
114        self.onlap_segments[:] = self.insert_onlap_surfaces()
115        # if self.cfg.include_channels:
116        #     floodplain_shale, channel_fill, shale_channel_drape, levee, crevasse = self.build_channel_cubes()
117        #     self.floodplain_shale = self.vertical_anti_alias_filter_simple(floodplain_shale)
118        #     self.channel_fill = self.vertical_anti_alias_filter_simple(channel_fill)
119        #     self.shale_channel_drape = self.vertical_anti_alias_filter_simple(shale_channel_drape)
120        #     self.levee = self.vertical_anti_alias_filter_simple(levee)
121        #     self.crevasse = self.vertical_anti_alias_filter_simple(crevasse)
122
123    def infilled_cube_shape(self):
124        cube_shape = (
125            self.cfg.cube_shape[0],
126            self.cfg.cube_shape[1],
127            (self.cfg.pad_samples + self.cfg.cube_shape[2]) * self.cfg.infill_factor,
128        )
129        return cube_shape
130
131    def build_meshgrid(self):
132        """
133            Build Meshgrid
134            --------------------------
135            Creates a meshgrid using the data in the Geomodel.
136
137            Parameters
138            ----------
139            None
140
141            Returns
142            -------
143            meshgrid : np.darray
144                A meshgrid of the data in the Geomodel.
145        """
146        return np.meshgrid(
147            range(self.cfg.cube_shape[0]),
148            range(self.cfg.cube_shape[1]),
149            sparse=False,
150            indexing="ij",
151        )
152
153    def create_geologic_age_3d_from_infilled_horizons(self, depth_maps, verbose=False):
154        """
155            Create geologic age 3d model from infilled horizons.
156            --------------------------
157            Creates cube containing a geologic age model from a horizon stack.
158
159            - depth_maps have units like ft or ms
160            - geologic age has arbitrary units where 'age' is same as horizon index
161
162            Parameters
163            ----------
164            depth_maps : np.ndarray
165                THe depth maps to use to generate the geologic age model.
166            verbose : bool
167                The level of verbosity in the logs
168
169            Returns
170            -------
171            returns : age : np.ndarray
172                The geologic age model.
173        """
174        if self.cfg.verbose:
175            print("\nCreating Geologic Age volume from unfaulted depth maps")
176        cube_shape = self.infilled_cube_shape()
177        # ensure that first depth_map has zeros at all X,Y locations
178        if not np.all(depth_maps[:, :, 0] == 0.0):
179            depth_maps_temp = np.dstack((np.zeros(cube_shape[:2], "float"), depth_maps))
180        else:
181            depth_maps_temp = depth_maps.copy()
182
183        if verbose:
184            print("\n\n   ... inside create_geologic_age_3D_from_infilled_horizons ")
185            print(
186                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
187                    depth_maps[:, :, :-1].min(),
188                    depth_maps[:, :, :-1].mean(),
189                    depth_maps[:, :, :-1].max(),
190                    cube_shape,
191                )
192            )
193
194        # create geologic age cube
195        age_range = np.linspace(0.0, float(cube_shape[2] - 1), cube_shape[2])
196        age = np.zeros(cube_shape, "float")
197        for i in range(cube_shape[0]):
198            for j in range(cube_shape[1]):
199                index_max_geo_age = np.argmax(
200                    depth_maps_temp[i, j, :].clip(0.0, float(cube_shape[2] - 1))
201                )
202                age[i, j, :] = np.interp(
203                    age_range,
204                    depth_maps_temp[i, j, : int(index_max_geo_age)],
205                    np.arange(index_max_geo_age),
206                )
207
208        if self.cfg.verbose:
209            print(f"    ... age.shape = {age.shape}")
210            print(
211                f"    ... age min/mean/max = {age[:].min()}, {age[:].mean():.1f}, {age[:].max()}"
212            )
213            print(
214                "    ... finished create_geologic_age_3D_from_infilled_horizons ...\n"
215            )
216
217            from datagenerator.util import plot_xsection
218            from datagenerator.util import plot_voxels_not_in_regular_layers
219
220            plot_xsection(
221                age[:],
222                depth_maps,
223                cfg=self.cfg,
224                line_num=int(cube_shape[0] / 2),
225                title=f"relative geologic age interpolated from horizons\nprior to faulting\nInline:"
226                f" {int(cube_shape[0] / 2)}",
227                png_name="QC_plot__unfaulted_relative_geologic_age.png",
228            )
229
230            # Analyse voxels not in regular layers
231            plot_voxels_not_in_regular_layers(
232                volume=age[:],
233                threshold=0.0,
234                cfg=self.cfg,
235                title="Example Trav through 3D model\nhistogram of raw layer values",
236                png_name="QC_plot__histogram_raw_Layers_1.png",
237            )
238
239        return self.vertical_anti_alias_filter_simple(age)
240
241    def vertical_anti_alias_filter_simple(self, cube) -> np.ndarray:
242        """
243            Simple vertical anti alias filter
244            --------------------------
245            Applies a simple version of the vertical anti-alias filter.
246
247            Parameters
248            ----------
249            cube : np.ndarray
250                The cube to apply the filter to.
251            Returns
252            -------
253            filtered_cube : np.ndarray
254                The anti-aliased filtered cube.
255        """
256        if self.cfg.verbose:
257            print("Applying simple vertical anti-alias filter")
258        filtered_cube = cube[..., :: self.cfg.infill_factor]
259        return filtered_cube
260
261    def vertical_anti_alias_filter(self, cube):
262        """
263            Vertical anti alias filter
264            --------------------------
265            Applies a vertical anti-alias filter.
266
267            Parameters
268            ----------
269            cube : np.ndarray
270                The cube to apply the filter to.
271            Returns
272            -------
273            filtered_cube : np.ndarray
274                The anti-aliased filtered cube.
275        """
276        infill_factor_o2_odd = next_odd(int(self.cfg.infill_factor + 1) / 2)
277        max_filt = maximum_filter(cube, size=(1, 1, infill_factor_o2_odd + 2))
278        filtered_cube = self.vertical_anti_alias_filter_simple(max_filt)
279        return filtered_cube
280
281    def plot_channels(self, channel_segments):
282        """
283        Plot Channels
284        -------------
285
286        Plot the channel segments.
287
288        Parameters
289        ----------
290        channel_segments : np.ndarray
291            An array of the channel segments.
292
293        Returns
294        -------
295        None
296        """
297        from datagenerator.util import plot_voxels_not_in_regular_layers
298        from datagenerator.util import find_line_with_most_voxels
299        from datagenerator.util import import_matplotlib
300
301        plt = import_matplotlib()
302
303        # find inline with the most channel voxels
304        inline_index = find_line_with_most_voxels(
305            channel_segments, (1000, 2000), self.cfg
306        )
307        plt.close(1)
308        plt.figure(1, figsize=(10, 7))
309        plt.clf()
310        plt.title(
311            "Example Trav through 3D model\nLayers filled with Channels Facies Codes"
312        )
313        plt.imshow(
314            np.fliplr(
315                np.rot90(
316                    channel_segments[
317                        inline_index,
318                        :,
319                        : self.cfg.cube_shape[2] * self.cfg.infill_factor - 1,
320                    ],
321                    3,
322                )
323            ),
324            aspect="auto",
325            cmap="jet",
326        )
327        plt.colorbar()
328        plt.ylim((self.geologic_age.shape[-1], 0))
329        for i in range(0, self.depth_maps.shape[-1], 5):
330            plt.plot(
331                range(self.cfg.cube_shape[0]),
332                self.depth_maps[inline_index, :, i],
333                "k-",
334                lw=0.3,
335            )
336        image_path = os.path.join(
337            self.cfg.work_subfolder, "QC_plot__LayersFilledWith_ChannelsFacies.png"
338        )
339        plt.savefig(image_path, format="png")
340        plt.close()
341
342        # analyze voxel values not in regular layers
343        title = "Example Trav through 3D model\nhistogram of voxels related to channel facies / before faulting"
344        pname = "QC_plot__Channels__histogram_beforeFaulting.png"
345        plot_voxels_not_in_regular_layers(
346            channel_segments, 200.0, title, pname, self.cfg
347        )
348
349        title = "Example Trav through 3D model\nhistogram of raw layer values\nShould NOT see channel facies"
350        pname = "QC_plot__histogram_raw_Layers_2.png"
351        plot_voxels_not_in_regular_layers(
352            self.geologic_age, 0.0, title, pname, self.cfg
353        )
354
355    def insert_onlap_surfaces(self):
356        """
357        Insert onlap surfaces
358        -------------
359
360        Insert onlap surfaces into the geomodel.
361
362        Parameters
363        ----------
364        None
365
366        Returns
367        -------
368        onlap_segments : np.ndarray
369            An array of the onlap segments.
370        """
371        if self.cfg.verbose:
372            print(
373                "\n\n ... create 3D (pre-faulting) labels for tilting episodes"
374                "\n  ... reminder: tilting events were added at horizons {}\n".format(
375                    self.onlap_horizon_list
376                )
377            )
378
379        # get work_cube_shape from cfg parameters
380        cube_shape = self.infilled_cube_shape()
381
382        # make sure
383        if self.cfg.verbose:
384            print("\n\n   ... inside insertOnlap3Dsurface_prefault ")
385            print(
386                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
387                    self.depth_maps[:].min(),
388                    self.depth_maps[:].mean(),
389                    self.depth_maps[:].max(),
390                    cube_shape,
391                )
392            )
393
394        # create 3D cube to hold segmentation results
395        onlap_segments = np.zeros(cube_shape, "float32")
396
397        # create grids with grid indices
398        ii, jj = np.meshgrid(
399            range(cube_shape[0]), range(cube_shape[1]), sparse=False, indexing="ij"
400        )
401
402        # loop through horizons in 'onlaps_horizon_list'
403        if len(self.onlap_horizon_list) == 0:
404            return self.vertical_anti_alias_filter_simple(onlap_segments)
405
406        sorted_onlaps_horizon_list = self.onlap_horizon_list
407        sorted_onlaps_horizon_list.sort()
408
409        voxel_count = 0
410        for ihorizon in sorted_onlaps_horizon_list:
411
412            # insert onlap label in voxels around horizon
413            shallow_map = self.depth_maps[:, :, ihorizon].copy()
414            shallow_depth_map_integer = (shallow_map + 0.5).astype("int")
415
416            for k in range(
417                -int(self.cfg.infill_factor * 1.5),
418                int(self.cfg.infill_factor * 1.5) + 1,
419            ):
420
421                sublayer_ii = ii
422                sublayer_jj = jj
423                sublayer_depth_map = k + shallow_depth_map_integer
424
425                sublayer_depth_map = np.clip(sublayer_depth_map, 0, cube_shape[2] - 1)
426                sublayer = onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map]
427                voxel_count += sublayer.flatten().shape[0]
428
429                del sublayer
430
431                onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map] += 1.0
432                if self.cfg.verbose and ihorizon % 1 == 0:
433                    print(
434                        "\t... k: {}, voxel_count: {}, sublayer_current_depth_map.mean: {:.2f} ".format(
435                            k, voxel_count, sublayer_depth_map.mean()
436                        )
437                    )
438
439        # reset invalid values
440        onlap_segments[onlap_segments < 0.0] = 0.0
441
442        non_zero_pixels = onlap_segments[onlap_segments != 0.0].shape[0]
443        pct_non_zero = float(non_zero_pixels) / (
444            cube_shape[0] * cube_shape[1] * cube_shape[2]
445        )
446        if self.cfg.verbose:
447            print(
448                "\t...onlap_segments min: {}, mean: {}, max: {}, % non-zero: {} ".format(
449                    onlap_segments.min(),
450                    onlap_segments.mean(),
451                    onlap_segments.max(),
452                    pct_non_zero,
453                )
454            )
455
456        if self.cfg.verbose:
457            print("    ... finished putting onlap surface in onlap_segments ...\n")
458
459        # Print non zero voxel count details
460        if self.cfg.verbose:
461            voxel_count_non_zero = onlap_segments[onlap_segments != 0].shape[0]
462            print(
463                f"\n   ...(pre-faulting) onlap segments created. min: {onlap_segments.min()},"
464                f" mean: {onlap_segments.mean():.5f}, max: {onlap_segments.max()},"
465                f" voxel count: {voxel_count_non_zero}\n"
466                f"\n   ...(pre-faulting) onlap segments created.shape = {onlap_segments.shape}"
467            )
468
469        return self.vertical_anti_alias_filter_simple(onlap_segments)
470
471    def write_cube_to_disk(self, data: np.ndarray, fname: str):
472        """
473        Writes cube to disk.
474        -------------
475        Writes a cube to disk.
476
477        Parameters
478        ----------
479        data : np.ndarray
480            The data to be written to disk.
481        fname : str
482            The name to be influded in the file name.
483
484        Returns
485        -------
486        None
487
488        It generates a `.npy` file on disk.
489        """
490        """Write 3D array to npy format."""
491        fname = os.path.join(
492            self.cfg.work_subfolder, f"{fname}_{self.cfg.date_stamp}.npy"
493        )
494        np.save(fname, data)
class Geomodel: 
  9class Geomodel:
 10    """
 11        Geomodel
 12        --------------------------
 13        The class of the Geomodel object.
 14
 15        This class contains all the items that make up the Geologic model.
 16
 17        Parameters
 18        ----------
 19        parameters : datagenerator.Parameters
 20            Parameter object storing all model parameters.
 21        depth_maps : np.ndarray
 22            A numpy array containing the depth maps.
 23        onlap_horizon_list : list
 24            A list of the onlap horizons.
 25        facies : np.ndarray
 26            The generated facies.
 27
 28        Returns
 29        -------
 30        None
 31    """
 32    def __init__(
 33        self,
 34        parameters: Parameters,
 35        depth_maps: np.ndarray,
 36        onlap_horizon_list: list,
 37        facies: np.ndarray
 38    ) -> None:
 39        """__init__
 40
 41        Initializer for the Geomodel class.
 42
 43        Parameters
 44        ----------
 45        parameters : datagenerator.Parameters
 46            Parameter object storing all model parameters.
 47        depth_maps : np.ndarray
 48            A numpy array containing the depth maps.
 49        onlap_horizon_list : list
 50            A list of the onlap horizons.
 51        facies : np.ndarray
 52            The generated facies.
 53        """
 54        self.cfg = parameters
 55        self.depth_maps = depth_maps
 56        self.onlap_horizon_list = onlap_horizon_list
 57        self.facies = facies
 58
 59        # Initialise geomodels
 60        cube_shape = (
 61            self.cfg.cube_shape[0],
 62            self.cfg.cube_shape[1],
 63            self.cfg.cube_shape[2] + self.cfg.pad_samples,
 64        )
 65        self.geologic_age = self.cfg.hdf_init("geologic_age_prefault", shape=cube_shape)
 66        self.onlap_segments = self.cfg.hdf_init(
 67            "onlap_segments_prefault", shape=cube_shape
 68        )
 69        self.faulted_lithology = self.cfg.hdf_init(
 70            "lithology_prefault", shape=cube_shape
 71        )
 72        self.geomodel_ng = self.cfg.hdf_init("net_to_gross_prefault", shape=cube_shape)
 73        self.faulted_depth = self.cfg.hdf_init("depth_prefault", shape=cube_shape)
 74        self.faulted_depth_randomised = self.cfg.hdf_init(
 75            "depth_randomised_prefault", shape=cube_shape
 76        )
 77
 78        # Channel volumes
 79        if self.cfg.include_channels:
 80            self.floodplain_shale = None
 81            self.channel_fill = None
 82            self.shale_channel_drape = None
 83            self.levee = None
 84            self.crevasse = None
 85            self.channel_segments = None
 86
 87    def build_unfaulted_geomodels(self):
 88        """
 89            Build unfaulted geomodels.
 90            --------------------------
 91            A method that builds unfaulted geomodels.
 92
 93            This method does the following:
 94
 95            * Build geologic age cube from depth horizons
 96            * Build onlap segmentation cube
 97            * If channels are turned on, use fluvsim fortran code to build:
 98            - floodplain shale cube
 99            - channel fill cube
100            - shale channel drape
101            - levee cube
102            - crevasse cube
103
104            Parameters
105            ----------
106            None
107
108            Returns
109            -------
110            None
111        """
112        self.geologic_age[:] = self.create_geologic_age_3d_from_infilled_horizons(
113            self.depth_maps
114        )
115        self.onlap_segments[:] = self.insert_onlap_surfaces()
116        # if self.cfg.include_channels:
117        #     floodplain_shale, channel_fill, shale_channel_drape, levee, crevasse = self.build_channel_cubes()
118        #     self.floodplain_shale = self.vertical_anti_alias_filter_simple(floodplain_shale)
119        #     self.channel_fill = self.vertical_anti_alias_filter_simple(channel_fill)
120        #     self.shale_channel_drape = self.vertical_anti_alias_filter_simple(shale_channel_drape)
121        #     self.levee = self.vertical_anti_alias_filter_simple(levee)
122        #     self.crevasse = self.vertical_anti_alias_filter_simple(crevasse)
123
124    def infilled_cube_shape(self):
125        cube_shape = (
126            self.cfg.cube_shape[0],
127            self.cfg.cube_shape[1],
128            (self.cfg.pad_samples + self.cfg.cube_shape[2]) * self.cfg.infill_factor,
129        )
130        return cube_shape
131
132    def build_meshgrid(self):
133        """
134            Build Meshgrid
135            --------------------------
136            Creates a meshgrid using the data in the Geomodel.
137
138            Parameters
139            ----------
140            None
141
142            Returns
143            -------
144            meshgrid : np.darray
145                A meshgrid of the data in the Geomodel.
146        """
147        return np.meshgrid(
148            range(self.cfg.cube_shape[0]),
149            range(self.cfg.cube_shape[1]),
150            sparse=False,
151            indexing="ij",
152        )
153
154    def create_geologic_age_3d_from_infilled_horizons(self, depth_maps, verbose=False):
155        """
156            Create geologic age 3d model from infilled horizons.
157            --------------------------
158            Creates cube containing a geologic age model from a horizon stack.
159
160            - depth_maps have units like ft or ms
161            - geologic age has arbitrary units where 'age' is same as horizon index
162
163            Parameters
164            ----------
165            depth_maps : np.ndarray
166                THe depth maps to use to generate the geologic age model.
167            verbose : bool
168                The level of verbosity in the logs
169
170            Returns
171            -------
172            returns : age : np.ndarray
173                The geologic age model.
174        """
175        if self.cfg.verbose:
176            print("\nCreating Geologic Age volume from unfaulted depth maps")
177        cube_shape = self.infilled_cube_shape()
178        # ensure that first depth_map has zeros at all X,Y locations
179        if not np.all(depth_maps[:, :, 0] == 0.0):
180            depth_maps_temp = np.dstack((np.zeros(cube_shape[:2], "float"), depth_maps))
181        else:
182            depth_maps_temp = depth_maps.copy()
183
184        if verbose:
185            print("\n\n   ... inside create_geologic_age_3D_from_infilled_horizons ")
186            print(
187                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
188                    depth_maps[:, :, :-1].min(),
189                    depth_maps[:, :, :-1].mean(),
190                    depth_maps[:, :, :-1].max(),
191                    cube_shape,
192                )
193            )
194
195        # create geologic age cube
196        age_range = np.linspace(0.0, float(cube_shape[2] - 1), cube_shape[2])
197        age = np.zeros(cube_shape, "float")
198        for i in range(cube_shape[0]):
199            for j in range(cube_shape[1]):
200                index_max_geo_age = np.argmax(
201                    depth_maps_temp[i, j, :].clip(0.0, float(cube_shape[2] - 1))
202                )
203                age[i, j, :] = np.interp(
204                    age_range,
205                    depth_maps_temp[i, j, : int(index_max_geo_age)],
206                    np.arange(index_max_geo_age),
207                )
208
209        if self.cfg.verbose:
210            print(f"    ... age.shape = {age.shape}")
211            print(
212                f"    ... age min/mean/max = {age[:].min()}, {age[:].mean():.1f}, {age[:].max()}"
213            )
214            print(
215                "    ... finished create_geologic_age_3D_from_infilled_horizons ...\n"
216            )
217
218            from datagenerator.util import plot_xsection
219            from datagenerator.util import plot_voxels_not_in_regular_layers
220
221            plot_xsection(
222                age[:],
223                depth_maps,
224                cfg=self.cfg,
225                line_num=int(cube_shape[0] / 2),
226                title=f"relative geologic age interpolated from horizons\nprior to faulting\nInline:"
227                f" {int(cube_shape[0] / 2)}",
228                png_name="QC_plot__unfaulted_relative_geologic_age.png",
229            )
230
231            # Analyse voxels not in regular layers
232            plot_voxels_not_in_regular_layers(
233                volume=age[:],
234                threshold=0.0,
235                cfg=self.cfg,
236                title="Example Trav through 3D model\nhistogram of raw layer values",
237                png_name="QC_plot__histogram_raw_Layers_1.png",
238            )
239
240        return self.vertical_anti_alias_filter_simple(age)
241
242    def vertical_anti_alias_filter_simple(self, cube) -> np.ndarray:
243        """
244            Simple vertical anti alias filter
245            --------------------------
246            Applies a simple version of the vertical anti-alias filter.
247
248            Parameters
249            ----------
250            cube : np.ndarray
251                The cube to apply the filter to.
252            Returns
253            -------
254            filtered_cube : np.ndarray
255                The anti-aliased filtered cube.
256        """
257        if self.cfg.verbose:
258            print("Applying simple vertical anti-alias filter")
259        filtered_cube = cube[..., :: self.cfg.infill_factor]
260        return filtered_cube
261
262    def vertical_anti_alias_filter(self, cube):
263        """
264            Vertical anti alias filter
265            --------------------------
266            Applies a vertical anti-alias filter.
267
268            Parameters
269            ----------
270            cube : np.ndarray
271                The cube to apply the filter to.
272            Returns
273            -------
274            filtered_cube : np.ndarray
275                The anti-aliased filtered cube.
276        """
277        infill_factor_o2_odd = next_odd(int(self.cfg.infill_factor + 1) / 2)
278        max_filt = maximum_filter(cube, size=(1, 1, infill_factor_o2_odd + 2))
279        filtered_cube = self.vertical_anti_alias_filter_simple(max_filt)
280        return filtered_cube
281
282    def plot_channels(self, channel_segments):
283        """
284        Plot Channels
285        -------------
286
287        Plot the channel segments.
288
289        Parameters
290        ----------
291        channel_segments : np.ndarray
292            An array of the channel segments.
293
294        Returns
295        -------
296        None
297        """
298        from datagenerator.util import plot_voxels_not_in_regular_layers
299        from datagenerator.util import find_line_with_most_voxels
300        from datagenerator.util import import_matplotlib
301
302        plt = import_matplotlib()
303
304        # find inline with the most channel voxels
305        inline_index = find_line_with_most_voxels(
306            channel_segments, (1000, 2000), self.cfg
307        )
308        plt.close(1)
309        plt.figure(1, figsize=(10, 7))
310        plt.clf()
311        plt.title(
312            "Example Trav through 3D model\nLayers filled with Channels Facies Codes"
313        )
314        plt.imshow(
315            np.fliplr(
316                np.rot90(
317                    channel_segments[
318                        inline_index,
319                        :,
320                        : self.cfg.cube_shape[2] * self.cfg.infill_factor - 1,
321                    ],
322                    3,
323                )
324            ),
325            aspect="auto",
326            cmap="jet",
327        )
328        plt.colorbar()
329        plt.ylim((self.geologic_age.shape[-1], 0))
330        for i in range(0, self.depth_maps.shape[-1], 5):
331            plt.plot(
332                range(self.cfg.cube_shape[0]),
333                self.depth_maps[inline_index, :, i],
334                "k-",
335                lw=0.3,
336            )
337        image_path = os.path.join(
338            self.cfg.work_subfolder, "QC_plot__LayersFilledWith_ChannelsFacies.png"
339        )
340        plt.savefig(image_path, format="png")
341        plt.close()
342
343        # analyze voxel values not in regular layers
344        title = "Example Trav through 3D model\nhistogram of voxels related to channel facies / before faulting"
345        pname = "QC_plot__Channels__histogram_beforeFaulting.png"
346        plot_voxels_not_in_regular_layers(
347            channel_segments, 200.0, title, pname, self.cfg
348        )
349
350        title = "Example Trav through 3D model\nhistogram of raw layer values\nShould NOT see channel facies"
351        pname = "QC_plot__histogram_raw_Layers_2.png"
352        plot_voxels_not_in_regular_layers(
353            self.geologic_age, 0.0, title, pname, self.cfg
354        )
355
356    def insert_onlap_surfaces(self):
357        """
358        Insert onlap surfaces
359        -------------
360
361        Insert onlap surfaces into the geomodel.
362
363        Parameters
364        ----------
365        None
366
367        Returns
368        -------
369        onlap_segments : np.ndarray
370            An array of the onlap segments.
371        """
372        if self.cfg.verbose:
373            print(
374                "\n\n ... create 3D (pre-faulting) labels for tilting episodes"
375                "\n  ... reminder: tilting events were added at horizons {}\n".format(
376                    self.onlap_horizon_list
377                )
378            )
379
380        # get work_cube_shape from cfg parameters
381        cube_shape = self.infilled_cube_shape()
382
383        # make sure
384        if self.cfg.verbose:
385            print("\n\n   ... inside insertOnlap3Dsurface_prefault ")
386            print(
387                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
388                    self.depth_maps[:].min(),
389                    self.depth_maps[:].mean(),
390                    self.depth_maps[:].max(),
391                    cube_shape,
392                )
393            )
394
395        # create 3D cube to hold segmentation results
396        onlap_segments = np.zeros(cube_shape, "float32")
397
398        # create grids with grid indices
399        ii, jj = np.meshgrid(
400            range(cube_shape[0]), range(cube_shape[1]), sparse=False, indexing="ij"
401        )
402
403        # loop through horizons in 'onlaps_horizon_list'
404        if len(self.onlap_horizon_list) == 0:
405            return self.vertical_anti_alias_filter_simple(onlap_segments)
406
407        sorted_onlaps_horizon_list = self.onlap_horizon_list
408        sorted_onlaps_horizon_list.sort()
409
410        voxel_count = 0
411        for ihorizon in sorted_onlaps_horizon_list:
412
413            # insert onlap label in voxels around horizon
414            shallow_map = self.depth_maps[:, :, ihorizon].copy()
415            shallow_depth_map_integer = (shallow_map + 0.5).astype("int")
416
417            for k in range(
418                -int(self.cfg.infill_factor * 1.5),
419                int(self.cfg.infill_factor * 1.5) + 1,
420            ):
421
422                sublayer_ii = ii
423                sublayer_jj = jj
424                sublayer_depth_map = k + shallow_depth_map_integer
425
426                sublayer_depth_map = np.clip(sublayer_depth_map, 0, cube_shape[2] - 1)
427                sublayer = onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map]
428                voxel_count += sublayer.flatten().shape[0]
429
430                del sublayer
431
432                onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map] += 1.0
433                if self.cfg.verbose and ihorizon % 1 == 0:
434                    print(
435                        "\t... k: {}, voxel_count: {}, sublayer_current_depth_map.mean: {:.2f} ".format(
436                            k, voxel_count, sublayer_depth_map.mean()
437                        )
438                    )
439
440        # reset invalid values
441        onlap_segments[onlap_segments < 0.0] = 0.0
442
443        non_zero_pixels = onlap_segments[onlap_segments != 0.0].shape[0]
444        pct_non_zero = float(non_zero_pixels) / (
445            cube_shape[0] * cube_shape[1] * cube_shape[2]
446        )
447        if self.cfg.verbose:
448            print(
449                "\t...onlap_segments min: {}, mean: {}, max: {}, % non-zero: {} ".format(
450                    onlap_segments.min(),
451                    onlap_segments.mean(),
452                    onlap_segments.max(),
453                    pct_non_zero,
454                )
455            )
456
457        if self.cfg.verbose:
458            print("    ... finished putting onlap surface in onlap_segments ...\n")
459
460        # Print non zero voxel count details
461        if self.cfg.verbose:
462            voxel_count_non_zero = onlap_segments[onlap_segments != 0].shape[0]
463            print(
464                f"\n   ...(pre-faulting) onlap segments created. min: {onlap_segments.min()},"
465                f" mean: {onlap_segments.mean():.5f}, max: {onlap_segments.max()},"
466                f" voxel count: {voxel_count_non_zero}\n"
467                f"\n   ...(pre-faulting) onlap segments created.shape = {onlap_segments.shape}"
468            )
469
470        return self.vertical_anti_alias_filter_simple(onlap_segments)
471
472    def write_cube_to_disk(self, data: np.ndarray, fname: str):
473        """
474        Writes cube to disk.
475        -------------
476        Writes a cube to disk.
477
478        Parameters
479        ----------
480        data : np.ndarray
481            The data to be written to disk.
482        fname : str
483            The name to be influded in the file name.
484
485        Returns
486        -------
487        None
488
489        It generates a `.npy` file on disk.
490        """
491        """Write 3D array to npy format."""
492        fname = os.path.join(
493            self.cfg.work_subfolder, f"{fname}_{self.cfg.date_stamp}.npy"
494        )
495        np.save(fname, data)
Geomodel

The class of the Geomodel object.

This class contains all the items that make up the Geologic model.

Parameters
  • parameters (datagenerator.Parameters): Parameter object storing all model parameters.
  • depth_maps (np.ndarray): A numpy array containing the depth maps.
  • onlap_horizon_list (list): A list of the onlap horizons.
  • facies (np.ndarray): The generated facies.
Returns
  • None
Geomodel( parameters: datagenerator.Parameters.Parameters, depth_maps: numpy.ndarray, onlap_horizon_list: list, facies: numpy.ndarray) 
32    def __init__(
33        self,
34        parameters: Parameters,
35        depth_maps: np.ndarray,
36        onlap_horizon_list: list,
37        facies: np.ndarray
38    ) -> None:
39        """__init__
40
41        Initializer for the Geomodel class.
42
43        Parameters
44        ----------
45        parameters : datagenerator.Parameters
46            Parameter object storing all model parameters.
47        depth_maps : np.ndarray
48            A numpy array containing the depth maps.
49        onlap_horizon_list : list
50            A list of the onlap horizons.
51        facies : np.ndarray
52            The generated facies.
53        """
54        self.cfg = parameters
55        self.depth_maps = depth_maps
56        self.onlap_horizon_list = onlap_horizon_list
57        self.facies = facies
58
59        # Initialise geomodels
60        cube_shape = (
61            self.cfg.cube_shape[0],
62            self.cfg.cube_shape[1],
63            self.cfg.cube_shape[2] + self.cfg.pad_samples,
64        )
65        self.geologic_age = self.cfg.hdf_init("geologic_age_prefault", shape=cube_shape)
66        self.onlap_segments = self.cfg.hdf_init(
67            "onlap_segments_prefault", shape=cube_shape
68        )
69        self.faulted_lithology = self.cfg.hdf_init(
70            "lithology_prefault", shape=cube_shape
71        )
72        self.geomodel_ng = self.cfg.hdf_init("net_to_gross_prefault", shape=cube_shape)
73        self.faulted_depth = self.cfg.hdf_init("depth_prefault", shape=cube_shape)
74        self.faulted_depth_randomised = self.cfg.hdf_init(
75            "depth_randomised_prefault", shape=cube_shape
76        )
77
78        # Channel volumes
79        if self.cfg.include_channels:
80            self.floodplain_shale = None
81            self.channel_fill = None
82            self.shale_channel_drape = None
83            self.levee = None
84            self.crevasse = None
85            self.channel_segments = None

__init__

Initializer for the Geomodel class.

Parameters
  • parameters (datagenerator.Parameters): Parameter object storing all model parameters.
  • depth_maps (np.ndarray): A numpy array containing the depth maps.
  • onlap_horizon_list (list): A list of the onlap horizons.
  • facies (np.ndarray): The generated facies.
def build_unfaulted_geomodels(self): 
 87    def build_unfaulted_geomodels(self):
 88        """
 89            Build unfaulted geomodels.
 90            --------------------------
 91            A method that builds unfaulted geomodels.
 92
 93            This method does the following:
 94
 95            * Build geologic age cube from depth horizons
 96            * Build onlap segmentation cube
 97            * If channels are turned on, use fluvsim fortran code to build:
 98            - floodplain shale cube
 99            - channel fill cube
100            - shale channel drape
101            - levee cube
102            - crevasse cube
103
104            Parameters
105            ----------
106            None
107
108            Returns
109            -------
110            None
111        """
112        self.geologic_age[:] = self.create_geologic_age_3d_from_infilled_horizons(
113            self.depth_maps
114        )
115        self.onlap_segments[:] = self.insert_onlap_surfaces()
116        # if self.cfg.include_channels:
117        #     floodplain_shale, channel_fill, shale_channel_drape, levee, crevasse = self.build_channel_cubes()
118        #     self.floodplain_shale = self.vertical_anti_alias_filter_simple(floodplain_shale)
119        #     self.channel_fill = self.vertical_anti_alias_filter_simple(channel_fill)
120        #     self.shale_channel_drape = self.vertical_anti_alias_filter_simple(shale_channel_drape)
121        #     self.levee = self.vertical_anti_alias_filter_simple(levee)
122        #     self.crevasse = self.vertical_anti_alias_filter_simple(crevasse)

Build unfaulted geomodels.

A method that builds unfaulted geomodels.

This method does the following:

  • Build geologic age cube from depth horizons
  • Build onlap segmentation cube
  • If channels are turned on, use fluvsim fortran code to build:
  • floodplain shale cube
  • channel fill cube
  • shale channel drape
  • levee cube
  • crevasse cube
Parameters
  • None
Returns
  • None
def infilled_cube_shape(self): 
124    def infilled_cube_shape(self):
125        cube_shape = (
126            self.cfg.cube_shape[0],
127            self.cfg.cube_shape[1],
128            (self.cfg.pad_samples + self.cfg.cube_shape[2]) * self.cfg.infill_factor,
129        )
130        return cube_shape
def build_meshgrid(self): 
132    def build_meshgrid(self):
133        """
134            Build Meshgrid
135            --------------------------
136            Creates a meshgrid using the data in the Geomodel.
137
138            Parameters
139            ----------
140            None
141
142            Returns
143            -------
144            meshgrid : np.darray
145                A meshgrid of the data in the Geomodel.
146        """
147        return np.meshgrid(
148            range(self.cfg.cube_shape[0]),
149            range(self.cfg.cube_shape[1]),
150            sparse=False,
151            indexing="ij",
152        )
Build Meshgrid

Creates a meshgrid using the data in the Geomodel.

Parameters
  • None
Returns
  • meshgrid (np.darray): A meshgrid of the data in the Geomodel.
def create_geologic_age_3d_from_infilled_horizons(self, depth_maps, verbose=False): 
154    def create_geologic_age_3d_from_infilled_horizons(self, depth_maps, verbose=False):
155        """
156            Create geologic age 3d model from infilled horizons.
157            --------------------------
158            Creates cube containing a geologic age model from a horizon stack.
159
160            - depth_maps have units like ft or ms
161            - geologic age has arbitrary units where 'age' is same as horizon index
162
163            Parameters
164            ----------
165            depth_maps : np.ndarray
166                THe depth maps to use to generate the geologic age model.
167            verbose : bool
168                The level of verbosity in the logs
169
170            Returns
171            -------
172            returns : age : np.ndarray
173                The geologic age model.
174        """
175        if self.cfg.verbose:
176            print("\nCreating Geologic Age volume from unfaulted depth maps")
177        cube_shape = self.infilled_cube_shape()
178        # ensure that first depth_map has zeros at all X,Y locations
179        if not np.all(depth_maps[:, :, 0] == 0.0):
180            depth_maps_temp = np.dstack((np.zeros(cube_shape[:2], "float"), depth_maps))
181        else:
182            depth_maps_temp = depth_maps.copy()
183
184        if verbose:
185            print("\n\n   ... inside create_geologic_age_3D_from_infilled_horizons ")
186            print(
187                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
188                    depth_maps[:, :, :-1].min(),
189                    depth_maps[:, :, :-1].mean(),
190                    depth_maps[:, :, :-1].max(),
191                    cube_shape,
192                )
193            )
194
195        # create geologic age cube
196        age_range = np.linspace(0.0, float(cube_shape[2] - 1), cube_shape[2])
197        age = np.zeros(cube_shape, "float")
198        for i in range(cube_shape[0]):
199            for j in range(cube_shape[1]):
200                index_max_geo_age = np.argmax(
201                    depth_maps_temp[i, j, :].clip(0.0, float(cube_shape[2] - 1))
202                )
203                age[i, j, :] = np.interp(
204                    age_range,
205                    depth_maps_temp[i, j, : int(index_max_geo_age)],
206                    np.arange(index_max_geo_age),
207                )
208
209        if self.cfg.verbose:
210            print(f"    ... age.shape = {age.shape}")
211            print(
212                f"    ... age min/mean/max = {age[:].min()}, {age[:].mean():.1f}, {age[:].max()}"
213            )
214            print(
215                "    ... finished create_geologic_age_3D_from_infilled_horizons ...\n"
216            )
217
218            from datagenerator.util import plot_xsection
219            from datagenerator.util import plot_voxels_not_in_regular_layers
220
221            plot_xsection(
222                age[:],
223                depth_maps,
224                cfg=self.cfg,
225                line_num=int(cube_shape[0] / 2),
226                title=f"relative geologic age interpolated from horizons\nprior to faulting\nInline:"
227                f" {int(cube_shape[0] / 2)}",
228                png_name="QC_plot__unfaulted_relative_geologic_age.png",
229            )
230
231            # Analyse voxels not in regular layers
232            plot_voxels_not_in_regular_layers(
233                volume=age[:],
234                threshold=0.0,
235                cfg=self.cfg,
236                title="Example Trav through 3D model\nhistogram of raw layer values",
237                png_name="QC_plot__histogram_raw_Layers_1.png",
238            )
239
240        return self.vertical_anti_alias_filter_simple(age)

Create geologic age 3d model from infilled horizons.

Creates cube containing a geologic age model from a horizon stack.

  • depth_maps have units like ft or ms
  • geologic age has arbitrary units where 'age' is same as horizon index
Parameters
  • depth_maps (np.ndarray): THe depth maps to use to generate the geologic age model.
  • verbose (bool): The level of verbosity in the logs
Returns
  • returns : age (np.ndarray): The geologic age model.
def vertical_anti_alias_filter_simple(self, cube) -> numpy.ndarray: 
242    def vertical_anti_alias_filter_simple(self, cube) -> np.ndarray:
243        """
244            Simple vertical anti alias filter
245            --------------------------
246            Applies a simple version of the vertical anti-alias filter.
247
248            Parameters
249            ----------
250            cube : np.ndarray
251                The cube to apply the filter to.
252            Returns
253            -------
254            filtered_cube : np.ndarray
255                The anti-aliased filtered cube.
256        """
257        if self.cfg.verbose:
258            print("Applying simple vertical anti-alias filter")
259        filtered_cube = cube[..., :: self.cfg.infill_factor]
260        return filtered_cube
Simple vertical anti alias filter

Applies a simple version of the vertical anti-alias filter.

Parameters
  • cube (np.ndarray): The cube to apply the filter to.
Returns
  • filtered_cube (np.ndarray): The anti-aliased filtered cube.
def vertical_anti_alias_filter(self, cube): 
262    def vertical_anti_alias_filter(self, cube):
263        """
264            Vertical anti alias filter
265            --------------------------
266            Applies a vertical anti-alias filter.
267
268            Parameters
269            ----------
270            cube : np.ndarray
271                The cube to apply the filter to.
272            Returns
273            -------
274            filtered_cube : np.ndarray
275                The anti-aliased filtered cube.
276        """
277        infill_factor_o2_odd = next_odd(int(self.cfg.infill_factor + 1) / 2)
278        max_filt = maximum_filter(cube, size=(1, 1, infill_factor_o2_odd + 2))
279        filtered_cube = self.vertical_anti_alias_filter_simple(max_filt)
280        return filtered_cube
Vertical anti alias filter

Applies a vertical anti-alias filter.

Parameters
  • cube (np.ndarray): The cube to apply the filter to.
Returns
  • filtered_cube (np.ndarray): The anti-aliased filtered cube.
def plot_channels(self, channel_segments): 
282    def plot_channels(self, channel_segments):
283        """
284        Plot Channels
285        -------------
286
287        Plot the channel segments.
288
289        Parameters
290        ----------
291        channel_segments : np.ndarray
292            An array of the channel segments.
293
294        Returns
295        -------
296        None
297        """
298        from datagenerator.util import plot_voxels_not_in_regular_layers
299        from datagenerator.util import find_line_with_most_voxels
300        from datagenerator.util import import_matplotlib
301
302        plt = import_matplotlib()
303
304        # find inline with the most channel voxels
305        inline_index = find_line_with_most_voxels(
306            channel_segments, (1000, 2000), self.cfg
307        )
308        plt.close(1)
309        plt.figure(1, figsize=(10, 7))
310        plt.clf()
311        plt.title(
312            "Example Trav through 3D model\nLayers filled with Channels Facies Codes"
313        )
314        plt.imshow(
315            np.fliplr(
316                np.rot90(
317                    channel_segments[
318                        inline_index,
319                        :,
320                        : self.cfg.cube_shape[2] * self.cfg.infill_factor - 1,
321                    ],
322                    3,
323                )
324            ),
325            aspect="auto",
326            cmap="jet",
327        )
328        plt.colorbar()
329        plt.ylim((self.geologic_age.shape[-1], 0))
330        for i in range(0, self.depth_maps.shape[-1], 5):
331            plt.plot(
332                range(self.cfg.cube_shape[0]),
333                self.depth_maps[inline_index, :, i],
334                "k-",
335                lw=0.3,
336            )
337        image_path = os.path.join(
338            self.cfg.work_subfolder, "QC_plot__LayersFilledWith_ChannelsFacies.png"
339        )
340        plt.savefig(image_path, format="png")
341        plt.close()
342
343        # analyze voxel values not in regular layers
344        title = "Example Trav through 3D model\nhistogram of voxels related to channel facies / before faulting"
345        pname = "QC_plot__Channels__histogram_beforeFaulting.png"
346        plot_voxels_not_in_regular_layers(
347            channel_segments, 200.0, title, pname, self.cfg
348        )
349
350        title = "Example Trav through 3D model\nhistogram of raw layer values\nShould NOT see channel facies"
351        pname = "QC_plot__histogram_raw_Layers_2.png"
352        plot_voxels_not_in_regular_layers(
353            self.geologic_age, 0.0, title, pname, self.cfg
354        )
Plot Channels

Plot the channel segments.

Parameters
  • channel_segments (np.ndarray): An array of the channel segments.
Returns
  • None
def insert_onlap_surfaces(self): 
356    def insert_onlap_surfaces(self):
357        """
358        Insert onlap surfaces
359        -------------
360
361        Insert onlap surfaces into the geomodel.
362
363        Parameters
364        ----------
365        None
366
367        Returns
368        -------
369        onlap_segments : np.ndarray
370            An array of the onlap segments.
371        """
372        if self.cfg.verbose:
373            print(
374                "\n\n ... create 3D (pre-faulting) labels for tilting episodes"
375                "\n  ... reminder: tilting events were added at horizons {}\n".format(
376                    self.onlap_horizon_list
377                )
378            )
379
380        # get work_cube_shape from cfg parameters
381        cube_shape = self.infilled_cube_shape()
382
383        # make sure
384        if self.cfg.verbose:
385            print("\n\n   ... inside insertOnlap3Dsurface_prefault ")
386            print(
387                "    ... depth_maps min/mean/max, cube_shape = {} {} {} {}".format(
388                    self.depth_maps[:].min(),
389                    self.depth_maps[:].mean(),
390                    self.depth_maps[:].max(),
391                    cube_shape,
392                )
393            )
394
395        # create 3D cube to hold segmentation results
396        onlap_segments = np.zeros(cube_shape, "float32")
397
398        # create grids with grid indices
399        ii, jj = np.meshgrid(
400            range(cube_shape[0]), range(cube_shape[1]), sparse=False, indexing="ij"
401        )
402
403        # loop through horizons in 'onlaps_horizon_list'
404        if len(self.onlap_horizon_list) == 0:
405            return self.vertical_anti_alias_filter_simple(onlap_segments)
406
407        sorted_onlaps_horizon_list = self.onlap_horizon_list
408        sorted_onlaps_horizon_list.sort()
409
410        voxel_count = 0
411        for ihorizon in sorted_onlaps_horizon_list:
412
413            # insert onlap label in voxels around horizon
414            shallow_map = self.depth_maps[:, :, ihorizon].copy()
415            shallow_depth_map_integer = (shallow_map + 0.5).astype("int")
416
417            for k in range(
418                -int(self.cfg.infill_factor * 1.5),
419                int(self.cfg.infill_factor * 1.5) + 1,
420            ):
421
422                sublayer_ii = ii
423                sublayer_jj = jj
424                sublayer_depth_map = k + shallow_depth_map_integer
425
426                sublayer_depth_map = np.clip(sublayer_depth_map, 0, cube_shape[2] - 1)
427                sublayer = onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map]
428                voxel_count += sublayer.flatten().shape[0]
429
430                del sublayer
431
432                onlap_segments[sublayer_ii, sublayer_jj, sublayer_depth_map] += 1.0
433                if self.cfg.verbose and ihorizon % 1 == 0:
434                    print(
435                        "\t... k: {}, voxel_count: {}, sublayer_current_depth_map.mean: {:.2f} ".format(
436                            k, voxel_count, sublayer_depth_map.mean()
437                        )
438                    )
439
440        # reset invalid values
441        onlap_segments[onlap_segments < 0.0] = 0.0
442
443        non_zero_pixels = onlap_segments[onlap_segments != 0.0].shape[0]
444        pct_non_zero = float(non_zero_pixels) / (
445            cube_shape[0] * cube_shape[1] * cube_shape[2]
446        )
447        if self.cfg.verbose:
448            print(
449                "\t...onlap_segments min: {}, mean: {}, max: {}, % non-zero: {} ".format(
450                    onlap_segments.min(),
451                    onlap_segments.mean(),
452                    onlap_segments.max(),
453                    pct_non_zero,
454                )
455            )
456
457        if self.cfg.verbose:
458            print("    ... finished putting onlap surface in onlap_segments ...\n")
459
460        # Print non zero voxel count details
461        if self.cfg.verbose:
462            voxel_count_non_zero = onlap_segments[onlap_segments != 0].shape[0]
463            print(
464                f"\n   ...(pre-faulting) onlap segments created. min: {onlap_segments.min()},"
465                f" mean: {onlap_segments.mean():.5f}, max: {onlap_segments.max()},"
466                f" voxel count: {voxel_count_non_zero}\n"
467                f"\n   ...(pre-faulting) onlap segments created.shape = {onlap_segments.shape}"
468            )
469
470        return self.vertical_anti_alias_filter_simple(onlap_segments)
Insert onlap surfaces

Insert onlap surfaces into the geomodel.

Parameters
  • None
Returns
  • onlap_segments (np.ndarray): An array of the onlap segments.
def write_cube_to_disk(self, data: numpy.ndarray, fname: str): 
472    def write_cube_to_disk(self, data: np.ndarray, fname: str):
473        """
474        Writes cube to disk.
475        -------------
476        Writes a cube to disk.
477
478        Parameters
479        ----------
480        data : np.ndarray
481            The data to be written to disk.
482        fname : str
483            The name to be influded in the file name.
484
485        Returns
486        -------
487        None
488
489        It generates a `.npy` file on disk.
490        """
491        """Write 3D array to npy format."""
492        fname = os.path.join(
493            self.cfg.work_subfolder, f"{fname}_{self.cfg.date_stamp}.npy"
494        )
495        np.save(fname, data)

Writes cube to disk.

Writes a cube to disk.

Parameters
  • data (np.ndarray): The data to be written to disk.
  • fname (str): The name to be influded in the file name.
Returns
  • None
  • It generates a .npy file on disk.