datagenerator.Salt
1from math import cos, sin 2import numpy as np 3from datagenerator.util import is_it_in_hull 4 5 6class SaltModel: 7 """ 8 Salt Model Class 9 ---------------- 10 11 This class creates a 3D salt body and inserts it into the input cube. 12 """ 13 def __init__(self, parameters) -> None: 14 """ 15 Initialization function 16 ----------------------- 17 18 Initializes the SaltModel class. 19 20 21 Parameters 22 ---------- 23 parameters : datagenerator.Parameters 24 The parameters of the project. 25 26 Returns 27 ------- 28 None 29 """ 30 self.cfg = parameters 31 cube_shape = ( 32 self.cfg.cube_shape[0], 33 self.cfg.cube_shape[1], 34 self.cfg.cube_shape[2] + self.cfg.pad_samples, 35 ) 36 self.salt_segments = self.cfg.hdf_init("salt_segments", shape=cube_shape) 37 self.points = [] 38 39 def compute_salt_body_segmentation(self) -> None: 40 """ 41 Compute Salt Body Segmentation 42 ------------------------------ 43 44 Creates a 3D salt body. 45 46 Parameters 47 ---------- 48 None 49 50 Returns 51 ------- 52 None 53 """ 54 print("\n\n...compute salt segments cube...") 55 salt_radius = np.random.triangular( 56 self.cfg.cube_shape[0] / 6, 57 self.cfg.cube_shape[0] / 5, 58 self.cfg.cube_shape[0] / 4, 59 ) 60 shallowest_salt = self.cfg.h5file.root.ModelData.faulted_depth_maps[ 61 :, :, 1 62 ].copy() 63 # shallowest_salt -= pad_samples / infill_factor 64 # shallowest_salt -= pad_samples 65 # shallowest_salt /= infill_factor 66 shallowest_salt = shallowest_salt[np.abs(shallowest_salt) < 1.0e10] 67 shallowest_salt = np.percentile(shallowest_salt[~np.isnan(shallowest_salt)], 99) 68 shallowest_salt += np.random.uniform(150, 300) 69 # output_cube = np.zeros(self.cfg.cube_shape) 70 print( 71 f" ...top of salt at {shallowest_salt}, with cube having shape {self.salt_segments.shape}" 72 ) 73 self.salt_segments[:] = self.insertSalt3D( 74 shallowest_salt, 75 salt_radius=salt_radius, 76 ) 77 78 @staticmethod 79 def create_circular_pointcloud(cx, cy, cz, radius, verbose=False): 80 """ 81 Create a circular pointcloud 82 ---------------------------- 83 84 Creates a circle of points with random jitter in (x, y, z). 85 86 87 Parameters 88 ---------- 89 cx : float 90 X-coordinate for centre of points 91 cy : float 92 Y-coordinate for centre of points 93 cz : float 94 Z-coordinate for centre of points 95 radius : float 96 Radius of circle 97 verbose : bool, optional 98 Print points while building circle. Defaults to False. 99 100 Returns 101 ------- 102 None 103 """ 104 points = [] 105 for iangle in np.arange(0.0, 360.0, 10.0): 106 xy_max_jitter = 0.1 * radius 107 x = ( 108 cx 109 + radius * cos(iangle * np.pi / 180.0) 110 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 111 ) 112 y = ( 113 cy 114 + radius * sin(iangle * np.pi / 180.0) 115 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 116 ) 117 z = cz + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 118 points.append([x, y, z]) 119 if verbose: 120 if iangle == 0.0: 121 print("\n") 122 print(f" ...angle, points = {iangle}, ({x}, {y}, {z})") 123 124 return points 125 126 def salt_circle_points( 127 self, 128 centre_points: list, 129 radii: list, 130 cx: float, 131 cy: float, 132 ) -> None: 133 """ 134 Salt Circle Points 135 ------------------ 136 137 Creates clouds of points representing the top or base of a salt body. 138 139 Parameters 140 ---------- 141 centre_points : list 142 Z-coordinates for deep, mid, shallow and tip locations 143 radii : list 144 Radius constants for deep, mid and shallow circles 145 cx : float) 146 X-coordinate for centre of point cloud 147 cy : float 148 Y-coordinate for centre of point cloud 149 150 Returns 151 ------- 152 None 153 """ 154 c_deep, c_mid, c_shallow, c_tip = centre_points 155 r_deep, r_mid, r_shallow = radii 156 print( 157 f"\n ...Centre_tip, Centre_shallow, Centre_mid, Centre_deep = {c_tip:.2f}, {c_shallow:.2f}, {c_mid:.2f}, {c_deep:.2f}" 158 ) 159 print( 160 f"\n ...radius_deep, radius_mid, radius_shallow = {r_deep:.2f}, {r_mid:.2f}, {r_shallow:.2f}" 161 ) 162 163 self.points += self.create_circular_pointcloud(cx, cy, c_deep, r_deep) 164 self.points += self.create_circular_pointcloud(cx, cy, c_mid, r_mid) 165 self.points += self.create_circular_pointcloud(cx, cy, c_shallow, r_shallow) 166 167 # Add points for the tip 168 xy_max_jitter = 0.1 * r_shallow 169 x = cx + np.random.uniform(-xy_max_jitter, xy_max_jitter) 170 y = cy + np.random.uniform(-xy_max_jitter, xy_max_jitter) 171 z = c_tip + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 172 self.points.append([x, y, z]) 173 print(f"\n ...tip (x, y, z) = ({x:.2f}, {y:.2f}, {z:.2f})") 174 175 def insertSalt3D(self, top, salt_radius=100.0) -> np.ndarray: 176 """ 177 Insert 3D Salt 178 -------------- 179 180 Generates a random salt body and insert into input_cube. 181 182 generate a random salt body and insert into input_cube 183 salt is modeled using 2 (almost) vertically aligned spheres connected 184 by a convex hull. Points inside the hull are considered salt. 185 center (in x,y) of upper and lower spheres are randomly jittered 186 so salt body is not always predominantly vertical 187 salt body is built using a 'bag of points'. Each point is 188 randomly jittered in x,y,z 189 190 Parameters 191 ---------- 192 top : float 193 Z-coordinate for top of salt body 194 salt_radius : float, optional 195 Radius of salt body. Defaults to 100.0. 196 197 Returns 198 ------- 199 None 200 """ 201 # Points built to model the crest of a salt body, with 202 # gentle anitclinal structure 203 # Z-coordinates for center of top salt 204 C1_deep = (top + salt_radius) * 1.1 205 C1_mid = C1_deep - salt_radius * np.random.uniform(0.35, 0.45) 206 C1_shallow = C1_deep - salt_radius * np.random.uniform(0.87, 0.93) 207 C1_tip = C1_deep - salt_radius * np.random.uniform(0.98, 1.02) 208 top_centres = [C1_deep, C1_mid, C1_shallow, C1_tip] 209 210 # Radius constants for top of salt 211 R1_deep = salt_radius * np.random.uniform(0.9, 1.1) 212 R1_mid = R1_deep * np.random.uniform(0.37, 0.43) 213 R1_shallow = R1_deep * np.random.uniform(0.18, 0.25) 214 top_radii = [R1_deep, R1_mid, R1_shallow] 215 216 cube_shape = self.salt_segments.shape 217 218 # x,y coords for top of salt 219 center_x = cube_shape[0] / 2 + np.random.uniform( 220 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 221 ) 222 center_y = cube_shape[1] / 2 + np.random.uniform( 223 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 224 ) 225 self.salt_circle_points(top_centres, top_radii, center_x, center_y) 226 227 # Build points to model the base of a salt body 228 C2_deep = min( 229 max( 230 C1_deep + salt_radius * np.random.uniform(2.3, 12.5), 231 cube_shape[2] - 2.0 * (salt_radius * np.random.uniform(-1.3, 3.5)), 232 ), 233 cube_shape[2] + 1.0 * (salt_radius), 234 ) 235 C2_mid = C2_deep + salt_radius * np.random.uniform(0.35, 0.45) / 2.0 236 C2_shallow = C2_deep + salt_radius * np.random.uniform(0.87, 0.93) / 2.0 237 C2_tip = C2_deep + salt_radius * np.random.uniform(0.98, 1.02) / 2.0 238 base_centres = [C2_deep, C2_mid, C2_shallow, C2_tip] 239 240 # radius constants for lower circle 241 R2_deep = R1_shallow / 2.0 * np.random.uniform(0.9, 1.1) 242 R2_mid = R2_deep * np.random.uniform(0.37, 0.43) 243 R2_shallow = R2_deep * np.random.uniform(0.18, 0.25) 244 base_radii = [R2_deep, R2_mid, R2_shallow] 245 246 # x,y coords for lower circle 247 x_center = center_x * np.random.uniform(0.3, 1.7) 248 y_center = center_y * np.random.uniform(0.3, 1.7) 249 250 self.salt_circle_points(base_centres, base_radii, x_center, y_center) 251 252 # build convex hull around all points for salt 253 254 # create list of indices of entire cube 255 xx, yy, zz = np.meshgrid( 256 list(range(cube_shape[0])), 257 list(range(cube_shape[1])), 258 list(range(cube_shape[2])), 259 sparse=False, 260 indexing="ij", 261 ) 262 xyz = ( 263 np.vstack((xx.flatten(), yy.flatten(), zz.flatten())) 264 .swapaxes(0, 1) 265 .astype("float") 266 ) 267 268 # check points in cube for inside or outside hull 269 in_hull = is_it_in_hull(self.points, xyz) 270 in_hull = in_hull.reshape(cube_shape) 271 272 salt_segments = np.zeros(cube_shape, "int16") 273 salt_segments[in_hull == True] = 1.0 274 275 return salt_segments 276 277 def update_depth_maps_with_salt_segments(self) -> None: 278 """ 279 Update Depth maps and DepthMapsGaps with salt segments 280 ------------------------------------------------------ 281 282 Updates the depth maps with salt segments 283 284 Parameters 285 ---------- 286 None 287 288 Returns 289 ------- 290 None 291 """ 292 ii, jj = np.meshgrid( 293 range(self.cfg.cube_shape[0]), 294 range(self.cfg.cube_shape[1]), 295 sparse=False, 296 indexing="ij", 297 ) 298 299 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 300 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 301 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 302 303 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 304 305 for ihorizon in range(0, depth_maps.shape[2] - 1): 306 print(" ...inserting salt in horizons...") 307 print( 308 " ...depth_maps min/mean/max = ", 309 depth_maps[:, :, ihorizon].min(), 310 depth_maps[:, :, ihorizon].mean(), 311 depth_maps[:, :, ihorizon].max(), 312 ) 313 depth_map = depth_maps[:, :, ihorizon].copy() 314 depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 315 316 # remove horizon picks inside salt (depth_maps_gaps) 317 depth_maps_gaps_horizon = depth_map_gaps.copy() 318 depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 319 320 _faulted_depth_map = depth_map.copy() 321 faulted_depth_map_indices = _faulted_depth_map.astype("int") 322 faulted_depth_map_indices = np.clip( 323 faulted_depth_map_indices, 324 0, 325 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 326 ) 327 328 _label = salt_segments[ii, jj, faulted_depth_map_indices] 329 depth_maps_gaps_horizon[ 330 depth_maps_gaps_horizon == 0.0 331 ] = np.nan # reset nan's 332 depth_maps_gaps_horizon[_label > 0] = np.nan 333 depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 334 try: 335 number_muted_points = depth_maps_gaps_horizon[_label > 0].shape[0] 336 if number_muted_points > 0: 337 print( 338 " ... horizon {} will have {} points muted inside salt".format( 339 ihorizon, number_muted_points 340 ) 341 ) 342 except: 343 pass 344 345 self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt 346 347 def update_depth_maps_with_salt_segments_drag(self): 348 """ 349 Update depth maps with salt segments frag 350 ----------------------------------------- 351 352 Updates depth maps with salt segments and drag. 353 354 Parameters 355 ---------- 356 None 357 358 Returns 359 ------- 360 depth_maps_salt : np.ndarray 361 The depth maps with salt. 362 depth_maps_gaps_salt : np.ndarray 363 The depth maps gaps with salt. 364 """ 365 from scipy import ndimage 366 367 ii, jj = np.meshgrid( 368 range(self.cfg.cube_shape[0]), 369 range(self.cfg.cube_shape[1]), 370 sparse=False, 371 indexing="ij", 372 ) 373 374 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 375 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 376 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 377 378 if self.cfg.model_qc_volumes: 379 np.save(f"{self.cfg.work_subfolder}/depth_maps_presalt.npy", depth_maps) 380 np.save( 381 f"{self.cfg.work_subfolder}/depth_maps_gaps_presalt.npy", 382 depth_maps_gaps, 383 ) 384 np.save(f"{self.cfg.work_subfolder}/salt_segments.npy", salt_segments) 385 386 depth_maps_salt = np.zeros_like(depth_maps_gaps) 387 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 388 389 relative_salt_depth = 0 390 391 for ihorizon in range(0, depth_maps.shape[2]): 392 print(" ...inserting salt in horizons...") 393 print( 394 " ...depth_maps min/mean/max = ", 395 depth_maps[:, :, ihorizon].min(), 396 depth_maps[:, :, ihorizon].mean(), 397 depth_maps[:, :, ihorizon].max(), 398 ) 399 depth_map = depth_maps[:, :, ihorizon].copy() 400 # depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 401 402 # remove horizon picks inside salt (depth_maps_gaps) 403 # depth_maps_gaps_horizon = depth_map_gaps.copy() 404 # depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 405 406 _faulted_depth_map = depth_map.copy() 407 faulted_depth_map_indices = _faulted_depth_map.astype("int") 408 faulted_depth_map_indices = np.clip( 409 faulted_depth_map_indices, 410 0, 411 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 412 ) 413 414 _label = salt_segments[ii, jj, faulted_depth_map_indices] 415 if np.max(_label) > 0: 416 relative_salt_depth += 1 417 418 # depth_maps_gaps_horizon[ 419 # depth_maps_gaps_horizon == 0.0 420 # ] = np.nan # reset nan's 421 422 # Apply vertical shift to the horizons where they are inside the salt body, relative to top of salt 423 # depth_maps_gaps_horizon[_label > 0] -= 2 * relative_salt_depth 424 # depth_maps_gaps_horizon = ndimage.gaussian_filter(depth_maps_gaps_horizon, 3) 425 426 # Do the same shift to the non-gapped horizon to build the facies volume 427 depth_map[_label > 0] -= 2 * relative_salt_depth 428 depth_map = ndimage.gaussian_filter(depth_map, 3) 429 430 # depth_maps_gaps_horizon[_label > 0] = np.nan 431 depth_maps_salt[..., ihorizon] = depth_map 432 433 # Re-apply blanking inside the salt to the gapped horizons 434 # _faulted_depth_map = depth_map.copy() 435 # faulted_depth_map_indices = _faulted_depth_map.astype("int") 436 # faulted_depth_map_indices = np.clip( 437 # faulted_depth_map_indices, 438 # 0, 439 # self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 440 # ) 441 # _label = salt_segments[ii, jj, faulted_depth_map_indices] 442 443 # depth_maps_gaps_horizon[_label > 0] = np.nan 444 # depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 445 446 try: 447 number_muted_points = depth_maps[_label > 0].shape[0] 448 if number_muted_points > 0: 449 print( 450 f" ... horizon {ihorizon} will have {number_muted_points} points muted inside salt" 451 ) 452 except: 453 pass 454 455 if self.cfg.model_qc_volumes: 456 np.save( 457 f"{self.cfg.work_subfolder}/depth_maps_salt_prepushdown.npy", 458 depth_maps_salt, 459 ) 460 461 depth_maps_salt = push_down_remove_negative_thickness(depth_maps_salt) 462 463 for ihorizon in range(0, depth_maps.shape[2] - 1): 464 depth_map_gaps_salt = depth_maps_salt[:, :, ihorizon].copy() 465 faulted_depth_map_indices = depth_map_gaps_salt.astype("int") 466 faulted_depth_map_indices = np.clip( 467 faulted_depth_map_indices, 468 0, 469 salt_segments.shape[2] - 1, 470 ) 471 _label = salt_segments[ii, jj, faulted_depth_map_indices] 472 depth_map_gaps_salt[_label > 0] = np.nan 473 depth_maps_gaps_salt[..., ihorizon] = depth_map_gaps_salt 474 475 depth_maps_gaps_salt[np.isnan(depth_maps_gaps)] = np.nan 476 477 # self.cfg.h5file.root.ModelData.faulted_depth_maps[:] = depth_maps 478 # self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt 479 480 if self.cfg.model_qc_volumes: 481 np.save(f"{self.cfg.work_subfolder}/depth_maps_salt.npy", depth_maps_salt) 482 np.save( 483 f"{self.cfg.work_subfolder}/depth_maps_gaps_salt.npy", 484 depth_maps_gaps_salt, 485 ) 486 np.save(f"{self.cfg.work_subfolder}/facies_label.npy", _label) 487 488 return depth_maps_salt, depth_maps_gaps_salt 489 490 491def push_down_remove_negative_thickness(depth_maps: np.ndarray) -> np.ndarray: 492 """ 493 Remove negative thicknesses 494 --------------------------- 495 496 Removes negative thicknesses. 497 498 Parameters 499 ---------- 500 depth_maps : np.ndarray 501 The depth maps. 502 503 Returns 504 ------- 505 depth_maps : np.ndarray 506 The depth maps with the negative thicknesses removed. 507 """ 508 for i in range(depth_maps.shape[-1] - 1, 1, -1): 509 layer_thickness = depth_maps[..., i] - depth_maps[..., i - 1] 510 if np.min(layer_thickness) < 0: 511 np.clip(layer_thickness, 0, a_max=None, out=layer_thickness) 512 depth_maps[..., i - 1] = depth_maps[..., i] - layer_thickness 513 514 return depth_maps
class
SaltModel:
7class SaltModel: 8 """ 9 Salt Model Class 10 ---------------- 11 12 This class creates a 3D salt body and inserts it into the input cube. 13 """ 14 def __init__(self, parameters) -> None: 15 """ 16 Initialization function 17 ----------------------- 18 19 Initializes the SaltModel class. 20 21 22 Parameters 23 ---------- 24 parameters : datagenerator.Parameters 25 The parameters of the project. 26 27 Returns 28 ------- 29 None 30 """ 31 self.cfg = parameters 32 cube_shape = ( 33 self.cfg.cube_shape[0], 34 self.cfg.cube_shape[1], 35 self.cfg.cube_shape[2] + self.cfg.pad_samples, 36 ) 37 self.salt_segments = self.cfg.hdf_init("salt_segments", shape=cube_shape) 38 self.points = [] 39 40 def compute_salt_body_segmentation(self) -> None: 41 """ 42 Compute Salt Body Segmentation 43 ------------------------------ 44 45 Creates a 3D salt body. 46 47 Parameters 48 ---------- 49 None 50 51 Returns 52 ------- 53 None 54 """ 55 print("\n\n...compute salt segments cube...") 56 salt_radius = np.random.triangular( 57 self.cfg.cube_shape[0] / 6, 58 self.cfg.cube_shape[0] / 5, 59 self.cfg.cube_shape[0] / 4, 60 ) 61 shallowest_salt = self.cfg.h5file.root.ModelData.faulted_depth_maps[ 62 :, :, 1 63 ].copy() 64 # shallowest_salt -= pad_samples / infill_factor 65 # shallowest_salt -= pad_samples 66 # shallowest_salt /= infill_factor 67 shallowest_salt = shallowest_salt[np.abs(shallowest_salt) < 1.0e10] 68 shallowest_salt = np.percentile(shallowest_salt[~np.isnan(shallowest_salt)], 99) 69 shallowest_salt += np.random.uniform(150, 300) 70 # output_cube = np.zeros(self.cfg.cube_shape) 71 print( 72 f" ...top of salt at {shallowest_salt}, with cube having shape {self.salt_segments.shape}" 73 ) 74 self.salt_segments[:] = self.insertSalt3D( 75 shallowest_salt, 76 salt_radius=salt_radius, 77 ) 78 79 @staticmethod 80 def create_circular_pointcloud(cx, cy, cz, radius, verbose=False): 81 """ 82 Create a circular pointcloud 83 ---------------------------- 84 85 Creates a circle of points with random jitter in (x, y, z). 86 87 88 Parameters 89 ---------- 90 cx : float 91 X-coordinate for centre of points 92 cy : float 93 Y-coordinate for centre of points 94 cz : float 95 Z-coordinate for centre of points 96 radius : float 97 Radius of circle 98 verbose : bool, optional 99 Print points while building circle. Defaults to False. 100 101 Returns 102 ------- 103 None 104 """ 105 points = [] 106 for iangle in np.arange(0.0, 360.0, 10.0): 107 xy_max_jitter = 0.1 * radius 108 x = ( 109 cx 110 + radius * cos(iangle * np.pi / 180.0) 111 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 112 ) 113 y = ( 114 cy 115 + radius * sin(iangle * np.pi / 180.0) 116 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 117 ) 118 z = cz + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 119 points.append([x, y, z]) 120 if verbose: 121 if iangle == 0.0: 122 print("\n") 123 print(f" ...angle, points = {iangle}, ({x}, {y}, {z})") 124 125 return points 126 127 def salt_circle_points( 128 self, 129 centre_points: list, 130 radii: list, 131 cx: float, 132 cy: float, 133 ) -> None: 134 """ 135 Salt Circle Points 136 ------------------ 137 138 Creates clouds of points representing the top or base of a salt body. 139 140 Parameters 141 ---------- 142 centre_points : list 143 Z-coordinates for deep, mid, shallow and tip locations 144 radii : list 145 Radius constants for deep, mid and shallow circles 146 cx : float) 147 X-coordinate for centre of point cloud 148 cy : float 149 Y-coordinate for centre of point cloud 150 151 Returns 152 ------- 153 None 154 """ 155 c_deep, c_mid, c_shallow, c_tip = centre_points 156 r_deep, r_mid, r_shallow = radii 157 print( 158 f"\n ...Centre_tip, Centre_shallow, Centre_mid, Centre_deep = {c_tip:.2f}, {c_shallow:.2f}, {c_mid:.2f}, {c_deep:.2f}" 159 ) 160 print( 161 f"\n ...radius_deep, radius_mid, radius_shallow = {r_deep:.2f}, {r_mid:.2f}, {r_shallow:.2f}" 162 ) 163 164 self.points += self.create_circular_pointcloud(cx, cy, c_deep, r_deep) 165 self.points += self.create_circular_pointcloud(cx, cy, c_mid, r_mid) 166 self.points += self.create_circular_pointcloud(cx, cy, c_shallow, r_shallow) 167 168 # Add points for the tip 169 xy_max_jitter = 0.1 * r_shallow 170 x = cx + np.random.uniform(-xy_max_jitter, xy_max_jitter) 171 y = cy + np.random.uniform(-xy_max_jitter, xy_max_jitter) 172 z = c_tip + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 173 self.points.append([x, y, z]) 174 print(f"\n ...tip (x, y, z) = ({x:.2f}, {y:.2f}, {z:.2f})") 175 176 def insertSalt3D(self, top, salt_radius=100.0) -> np.ndarray: 177 """ 178 Insert 3D Salt 179 -------------- 180 181 Generates a random salt body and insert into input_cube. 182 183 generate a random salt body and insert into input_cube 184 salt is modeled using 2 (almost) vertically aligned spheres connected 185 by a convex hull. Points inside the hull are considered salt. 186 center (in x,y) of upper and lower spheres are randomly jittered 187 so salt body is not always predominantly vertical 188 salt body is built using a 'bag of points'. Each point is 189 randomly jittered in x,y,z 190 191 Parameters 192 ---------- 193 top : float 194 Z-coordinate for top of salt body 195 salt_radius : float, optional 196 Radius of salt body. Defaults to 100.0. 197 198 Returns 199 ------- 200 None 201 """ 202 # Points built to model the crest of a salt body, with 203 # gentle anitclinal structure 204 # Z-coordinates for center of top salt 205 C1_deep = (top + salt_radius) * 1.1 206 C1_mid = C1_deep - salt_radius * np.random.uniform(0.35, 0.45) 207 C1_shallow = C1_deep - salt_radius * np.random.uniform(0.87, 0.93) 208 C1_tip = C1_deep - salt_radius * np.random.uniform(0.98, 1.02) 209 top_centres = [C1_deep, C1_mid, C1_shallow, C1_tip] 210 211 # Radius constants for top of salt 212 R1_deep = salt_radius * np.random.uniform(0.9, 1.1) 213 R1_mid = R1_deep * np.random.uniform(0.37, 0.43) 214 R1_shallow = R1_deep * np.random.uniform(0.18, 0.25) 215 top_radii = [R1_deep, R1_mid, R1_shallow] 216 217 cube_shape = self.salt_segments.shape 218 219 # x,y coords for top of salt 220 center_x = cube_shape[0] / 2 + np.random.uniform( 221 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 222 ) 223 center_y = cube_shape[1] / 2 + np.random.uniform( 224 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 225 ) 226 self.salt_circle_points(top_centres, top_radii, center_x, center_y) 227 228 # Build points to model the base of a salt body 229 C2_deep = min( 230 max( 231 C1_deep + salt_radius * np.random.uniform(2.3, 12.5), 232 cube_shape[2] - 2.0 * (salt_radius * np.random.uniform(-1.3, 3.5)), 233 ), 234 cube_shape[2] + 1.0 * (salt_radius), 235 ) 236 C2_mid = C2_deep + salt_radius * np.random.uniform(0.35, 0.45) / 2.0 237 C2_shallow = C2_deep + salt_radius * np.random.uniform(0.87, 0.93) / 2.0 238 C2_tip = C2_deep + salt_radius * np.random.uniform(0.98, 1.02) / 2.0 239 base_centres = [C2_deep, C2_mid, C2_shallow, C2_tip] 240 241 # radius constants for lower circle 242 R2_deep = R1_shallow / 2.0 * np.random.uniform(0.9, 1.1) 243 R2_mid = R2_deep * np.random.uniform(0.37, 0.43) 244 R2_shallow = R2_deep * np.random.uniform(0.18, 0.25) 245 base_radii = [R2_deep, R2_mid, R2_shallow] 246 247 # x,y coords for lower circle 248 x_center = center_x * np.random.uniform(0.3, 1.7) 249 y_center = center_y * np.random.uniform(0.3, 1.7) 250 251 self.salt_circle_points(base_centres, base_radii, x_center, y_center) 252 253 # build convex hull around all points for salt 254 255 # create list of indices of entire cube 256 xx, yy, zz = np.meshgrid( 257 list(range(cube_shape[0])), 258 list(range(cube_shape[1])), 259 list(range(cube_shape[2])), 260 sparse=False, 261 indexing="ij", 262 ) 263 xyz = ( 264 np.vstack((xx.flatten(), yy.flatten(), zz.flatten())) 265 .swapaxes(0, 1) 266 .astype("float") 267 ) 268 269 # check points in cube for inside or outside hull 270 in_hull = is_it_in_hull(self.points, xyz) 271 in_hull = in_hull.reshape(cube_shape) 272 273 salt_segments = np.zeros(cube_shape, "int16") 274 salt_segments[in_hull == True] = 1.0 275 276 return salt_segments 277 278 def update_depth_maps_with_salt_segments(self) -> None: 279 """ 280 Update Depth maps and DepthMapsGaps with salt segments 281 ------------------------------------------------------ 282 283 Updates the depth maps with salt segments 284 285 Parameters 286 ---------- 287 None 288 289 Returns 290 ------- 291 None 292 """ 293 ii, jj = np.meshgrid( 294 range(self.cfg.cube_shape[0]), 295 range(self.cfg.cube_shape[1]), 296 sparse=False, 297 indexing="ij", 298 ) 299 300 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 301 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 302 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 303 304 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 305 306 for ihorizon in range(0, depth_maps.shape[2] - 1): 307 print(" ...inserting salt in horizons...") 308 print( 309 " ...depth_maps min/mean/max = ", 310 depth_maps[:, :, ihorizon].min(), 311 depth_maps[:, :, ihorizon].mean(), 312 depth_maps[:, :, ihorizon].max(), 313 ) 314 depth_map = depth_maps[:, :, ihorizon].copy() 315 depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 316 317 # remove horizon picks inside salt (depth_maps_gaps) 318 depth_maps_gaps_horizon = depth_map_gaps.copy() 319 depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 320 321 _faulted_depth_map = depth_map.copy() 322 faulted_depth_map_indices = _faulted_depth_map.astype("int") 323 faulted_depth_map_indices = np.clip( 324 faulted_depth_map_indices, 325 0, 326 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 327 ) 328 329 _label = salt_segments[ii, jj, faulted_depth_map_indices] 330 depth_maps_gaps_horizon[ 331 depth_maps_gaps_horizon == 0.0 332 ] = np.nan # reset nan's 333 depth_maps_gaps_horizon[_label > 0] = np.nan 334 depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 335 try: 336 number_muted_points = depth_maps_gaps_horizon[_label > 0].shape[0] 337 if number_muted_points > 0: 338 print( 339 " ... horizon {} will have {} points muted inside salt".format( 340 ihorizon, number_muted_points 341 ) 342 ) 343 except: 344 pass 345 346 self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt 347 348 def update_depth_maps_with_salt_segments_drag(self): 349 """ 350 Update depth maps with salt segments frag 351 ----------------------------------------- 352 353 Updates depth maps with salt segments and drag. 354 355 Parameters 356 ---------- 357 None 358 359 Returns 360 ------- 361 depth_maps_salt : np.ndarray 362 The depth maps with salt. 363 depth_maps_gaps_salt : np.ndarray 364 The depth maps gaps with salt. 365 """ 366 from scipy import ndimage 367 368 ii, jj = np.meshgrid( 369 range(self.cfg.cube_shape[0]), 370 range(self.cfg.cube_shape[1]), 371 sparse=False, 372 indexing="ij", 373 ) 374 375 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 376 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 377 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 378 379 if self.cfg.model_qc_volumes: 380 np.save(f"{self.cfg.work_subfolder}/depth_maps_presalt.npy", depth_maps) 381 np.save( 382 f"{self.cfg.work_subfolder}/depth_maps_gaps_presalt.npy", 383 depth_maps_gaps, 384 ) 385 np.save(f"{self.cfg.work_subfolder}/salt_segments.npy", salt_segments) 386 387 depth_maps_salt = np.zeros_like(depth_maps_gaps) 388 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 389 390 relative_salt_depth = 0 391 392 for ihorizon in range(0, depth_maps.shape[2]): 393 print(" ...inserting salt in horizons...") 394 print( 395 " ...depth_maps min/mean/max = ", 396 depth_maps[:, :, ihorizon].min(), 397 depth_maps[:, :, ihorizon].mean(), 398 depth_maps[:, :, ihorizon].max(), 399 ) 400 depth_map = depth_maps[:, :, ihorizon].copy() 401 # depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 402 403 # remove horizon picks inside salt (depth_maps_gaps) 404 # depth_maps_gaps_horizon = depth_map_gaps.copy() 405 # depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 406 407 _faulted_depth_map = depth_map.copy() 408 faulted_depth_map_indices = _faulted_depth_map.astype("int") 409 faulted_depth_map_indices = np.clip( 410 faulted_depth_map_indices, 411 0, 412 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 413 ) 414 415 _label = salt_segments[ii, jj, faulted_depth_map_indices] 416 if np.max(_label) > 0: 417 relative_salt_depth += 1 418 419 # depth_maps_gaps_horizon[ 420 # depth_maps_gaps_horizon == 0.0 421 # ] = np.nan # reset nan's 422 423 # Apply vertical shift to the horizons where they are inside the salt body, relative to top of salt 424 # depth_maps_gaps_horizon[_label > 0] -= 2 * relative_salt_depth 425 # depth_maps_gaps_horizon = ndimage.gaussian_filter(depth_maps_gaps_horizon, 3) 426 427 # Do the same shift to the non-gapped horizon to build the facies volume 428 depth_map[_label > 0] -= 2 * relative_salt_depth 429 depth_map = ndimage.gaussian_filter(depth_map, 3) 430 431 # depth_maps_gaps_horizon[_label > 0] = np.nan 432 depth_maps_salt[..., ihorizon] = depth_map 433 434 # Re-apply blanking inside the salt to the gapped horizons 435 # _faulted_depth_map = depth_map.copy() 436 # faulted_depth_map_indices = _faulted_depth_map.astype("int") 437 # faulted_depth_map_indices = np.clip( 438 # faulted_depth_map_indices, 439 # 0, 440 # self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 441 # ) 442 # _label = salt_segments[ii, jj, faulted_depth_map_indices] 443 444 # depth_maps_gaps_horizon[_label > 0] = np.nan 445 # depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 446 447 try: 448 number_muted_points = depth_maps[_label > 0].shape[0] 449 if number_muted_points > 0: 450 print( 451 f" ... horizon {ihorizon} will have {number_muted_points} points muted inside salt" 452 ) 453 except: 454 pass 455 456 if self.cfg.model_qc_volumes: 457 np.save( 458 f"{self.cfg.work_subfolder}/depth_maps_salt_prepushdown.npy", 459 depth_maps_salt, 460 ) 461 462 depth_maps_salt = push_down_remove_negative_thickness(depth_maps_salt) 463 464 for ihorizon in range(0, depth_maps.shape[2] - 1): 465 depth_map_gaps_salt = depth_maps_salt[:, :, ihorizon].copy() 466 faulted_depth_map_indices = depth_map_gaps_salt.astype("int") 467 faulted_depth_map_indices = np.clip( 468 faulted_depth_map_indices, 469 0, 470 salt_segments.shape[2] - 1, 471 ) 472 _label = salt_segments[ii, jj, faulted_depth_map_indices] 473 depth_map_gaps_salt[_label > 0] = np.nan 474 depth_maps_gaps_salt[..., ihorizon] = depth_map_gaps_salt 475 476 depth_maps_gaps_salt[np.isnan(depth_maps_gaps)] = np.nan 477 478 # self.cfg.h5file.root.ModelData.faulted_depth_maps[:] = depth_maps 479 # self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt 480 481 if self.cfg.model_qc_volumes: 482 np.save(f"{self.cfg.work_subfolder}/depth_maps_salt.npy", depth_maps_salt) 483 np.save( 484 f"{self.cfg.work_subfolder}/depth_maps_gaps_salt.npy", 485 depth_maps_gaps_salt, 486 ) 487 np.save(f"{self.cfg.work_subfolder}/facies_label.npy", _label) 488 489 return depth_maps_salt, depth_maps_gaps_salt
Salt Model Class
This class creates a 3D salt body and inserts it into the input cube.
SaltModel(parameters)
14 def __init__(self, parameters) -> None: 15 """ 16 Initialization function 17 ----------------------- 18 19 Initializes the SaltModel class. 20 21 22 Parameters 23 ---------- 24 parameters : datagenerator.Parameters 25 The parameters of the project. 26 27 Returns 28 ------- 29 None 30 """ 31 self.cfg = parameters 32 cube_shape = ( 33 self.cfg.cube_shape[0], 34 self.cfg.cube_shape[1], 35 self.cfg.cube_shape[2] + self.cfg.pad_samples, 36 ) 37 self.salt_segments = self.cfg.hdf_init("salt_segments", shape=cube_shape) 38 self.points = []
Initialization function
Initializes the SaltModel class.
Parameters
- parameters (datagenerator.Parameters): The parameters of the project.
Returns
- None
def
compute_salt_body_segmentation(self) -> None:
40 def compute_salt_body_segmentation(self) -> None: 41 """ 42 Compute Salt Body Segmentation 43 ------------------------------ 44 45 Creates a 3D salt body. 46 47 Parameters 48 ---------- 49 None 50 51 Returns 52 ------- 53 None 54 """ 55 print("\n\n...compute salt segments cube...") 56 salt_radius = np.random.triangular( 57 self.cfg.cube_shape[0] / 6, 58 self.cfg.cube_shape[0] / 5, 59 self.cfg.cube_shape[0] / 4, 60 ) 61 shallowest_salt = self.cfg.h5file.root.ModelData.faulted_depth_maps[ 62 :, :, 1 63 ].copy() 64 # shallowest_salt -= pad_samples / infill_factor 65 # shallowest_salt -= pad_samples 66 # shallowest_salt /= infill_factor 67 shallowest_salt = shallowest_salt[np.abs(shallowest_salt) < 1.0e10] 68 shallowest_salt = np.percentile(shallowest_salt[~np.isnan(shallowest_salt)], 99) 69 shallowest_salt += np.random.uniform(150, 300) 70 # output_cube = np.zeros(self.cfg.cube_shape) 71 print( 72 f" ...top of salt at {shallowest_salt}, with cube having shape {self.salt_segments.shape}" 73 ) 74 self.salt_segments[:] = self.insertSalt3D( 75 shallowest_salt, 76 salt_radius=salt_radius, 77 )
Compute Salt Body Segmentation
Creates a 3D salt body.
Parameters
- None
Returns
- None
@staticmethod
def
create_circular_pointcloud(cx, cy, cz, radius, verbose=False):
79 @staticmethod 80 def create_circular_pointcloud(cx, cy, cz, radius, verbose=False): 81 """ 82 Create a circular pointcloud 83 ---------------------------- 84 85 Creates a circle of points with random jitter in (x, y, z). 86 87 88 Parameters 89 ---------- 90 cx : float 91 X-coordinate for centre of points 92 cy : float 93 Y-coordinate for centre of points 94 cz : float 95 Z-coordinate for centre of points 96 radius : float 97 Radius of circle 98 verbose : bool, optional 99 Print points while building circle. Defaults to False. 100 101 Returns 102 ------- 103 None 104 """ 105 points = [] 106 for iangle in np.arange(0.0, 360.0, 10.0): 107 xy_max_jitter = 0.1 * radius 108 x = ( 109 cx 110 + radius * cos(iangle * np.pi / 180.0) 111 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 112 ) 113 y = ( 114 cy 115 + radius * sin(iangle * np.pi / 180.0) 116 + np.random.uniform(-xy_max_jitter, xy_max_jitter) 117 ) 118 z = cz + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 119 points.append([x, y, z]) 120 if verbose: 121 if iangle == 0.0: 122 print("\n") 123 print(f" ...angle, points = {iangle}, ({x}, {y}, {z})") 124 125 return points
Create a circular pointcloud
Creates a circle of points with random jitter in (x, y, z).
Parameters
- cx (float): X-coordinate for centre of points
- cy (float): Y-coordinate for centre of points
- cz (float): Z-coordinate for centre of points
- radius (float): Radius of circle
- verbose (bool, optional): Print points while building circle. Defaults to False.
Returns
- None
def
salt_circle_points(self, centre_points: list, radii: list, cx: float, cy: float) -> None:
127 def salt_circle_points( 128 self, 129 centre_points: list, 130 radii: list, 131 cx: float, 132 cy: float, 133 ) -> None: 134 """ 135 Salt Circle Points 136 ------------------ 137 138 Creates clouds of points representing the top or base of a salt body. 139 140 Parameters 141 ---------- 142 centre_points : list 143 Z-coordinates for deep, mid, shallow and tip locations 144 radii : list 145 Radius constants for deep, mid and shallow circles 146 cx : float) 147 X-coordinate for centre of point cloud 148 cy : float 149 Y-coordinate for centre of point cloud 150 151 Returns 152 ------- 153 None 154 """ 155 c_deep, c_mid, c_shallow, c_tip = centre_points 156 r_deep, r_mid, r_shallow = radii 157 print( 158 f"\n ...Centre_tip, Centre_shallow, Centre_mid, Centre_deep = {c_tip:.2f}, {c_shallow:.2f}, {c_mid:.2f}, {c_deep:.2f}" 159 ) 160 print( 161 f"\n ...radius_deep, radius_mid, radius_shallow = {r_deep:.2f}, {r_mid:.2f}, {r_shallow:.2f}" 162 ) 163 164 self.points += self.create_circular_pointcloud(cx, cy, c_deep, r_deep) 165 self.points += self.create_circular_pointcloud(cx, cy, c_mid, r_mid) 166 self.points += self.create_circular_pointcloud(cx, cy, c_shallow, r_shallow) 167 168 # Add points for the tip 169 xy_max_jitter = 0.1 * r_shallow 170 x = cx + np.random.uniform(-xy_max_jitter, xy_max_jitter) 171 y = cy + np.random.uniform(-xy_max_jitter, xy_max_jitter) 172 z = c_tip + np.random.uniform(-xy_max_jitter / 2.0, xy_max_jitter / 2.0) 173 self.points.append([x, y, z]) 174 print(f"\n ...tip (x, y, z) = ({x:.2f}, {y:.2f}, {z:.2f})")
Salt Circle Points
Creates clouds of points representing the top or base of a salt body.
Parameters
- centre_points (list): Z-coordinates for deep, mid, shallow and tip locations
- radii (list): Radius constants for deep, mid and shallow circles
- cx (float)): X-coordinate for centre of point cloud
- cy (float): Y-coordinate for centre of point cloud
Returns
- None
def
insertSalt3D(self, top, salt_radius=100.0) -> numpy.ndarray:
176 def insertSalt3D(self, top, salt_radius=100.0) -> np.ndarray: 177 """ 178 Insert 3D Salt 179 -------------- 180 181 Generates a random salt body and insert into input_cube. 182 183 generate a random salt body and insert into input_cube 184 salt is modeled using 2 (almost) vertically aligned spheres connected 185 by a convex hull. Points inside the hull are considered salt. 186 center (in x,y) of upper and lower spheres are randomly jittered 187 so salt body is not always predominantly vertical 188 salt body is built using a 'bag of points'. Each point is 189 randomly jittered in x,y,z 190 191 Parameters 192 ---------- 193 top : float 194 Z-coordinate for top of salt body 195 salt_radius : float, optional 196 Radius of salt body. Defaults to 100.0. 197 198 Returns 199 ------- 200 None 201 """ 202 # Points built to model the crest of a salt body, with 203 # gentle anitclinal structure 204 # Z-coordinates for center of top salt 205 C1_deep = (top + salt_radius) * 1.1 206 C1_mid = C1_deep - salt_radius * np.random.uniform(0.35, 0.45) 207 C1_shallow = C1_deep - salt_radius * np.random.uniform(0.87, 0.93) 208 C1_tip = C1_deep - salt_radius * np.random.uniform(0.98, 1.02) 209 top_centres = [C1_deep, C1_mid, C1_shallow, C1_tip] 210 211 # Radius constants for top of salt 212 R1_deep = salt_radius * np.random.uniform(0.9, 1.1) 213 R1_mid = R1_deep * np.random.uniform(0.37, 0.43) 214 R1_shallow = R1_deep * np.random.uniform(0.18, 0.25) 215 top_radii = [R1_deep, R1_mid, R1_shallow] 216 217 cube_shape = self.salt_segments.shape 218 219 # x,y coords for top of salt 220 center_x = cube_shape[0] / 2 + np.random.uniform( 221 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 222 ) 223 center_y = cube_shape[1] / 2 + np.random.uniform( 224 -cube_shape[0] * 0.4, cube_shape[0] * 0.4 225 ) 226 self.salt_circle_points(top_centres, top_radii, center_x, center_y) 227 228 # Build points to model the base of a salt body 229 C2_deep = min( 230 max( 231 C1_deep + salt_radius * np.random.uniform(2.3, 12.5), 232 cube_shape[2] - 2.0 * (salt_radius * np.random.uniform(-1.3, 3.5)), 233 ), 234 cube_shape[2] + 1.0 * (salt_radius), 235 ) 236 C2_mid = C2_deep + salt_radius * np.random.uniform(0.35, 0.45) / 2.0 237 C2_shallow = C2_deep + salt_radius * np.random.uniform(0.87, 0.93) / 2.0 238 C2_tip = C2_deep + salt_radius * np.random.uniform(0.98, 1.02) / 2.0 239 base_centres = [C2_deep, C2_mid, C2_shallow, C2_tip] 240 241 # radius constants for lower circle 242 R2_deep = R1_shallow / 2.0 * np.random.uniform(0.9, 1.1) 243 R2_mid = R2_deep * np.random.uniform(0.37, 0.43) 244 R2_shallow = R2_deep * np.random.uniform(0.18, 0.25) 245 base_radii = [R2_deep, R2_mid, R2_shallow] 246 247 # x,y coords for lower circle 248 x_center = center_x * np.random.uniform(0.3, 1.7) 249 y_center = center_y * np.random.uniform(0.3, 1.7) 250 251 self.salt_circle_points(base_centres, base_radii, x_center, y_center) 252 253 # build convex hull around all points for salt 254 255 # create list of indices of entire cube 256 xx, yy, zz = np.meshgrid( 257 list(range(cube_shape[0])), 258 list(range(cube_shape[1])), 259 list(range(cube_shape[2])), 260 sparse=False, 261 indexing="ij", 262 ) 263 xyz = ( 264 np.vstack((xx.flatten(), yy.flatten(), zz.flatten())) 265 .swapaxes(0, 1) 266 .astype("float") 267 ) 268 269 # check points in cube for inside or outside hull 270 in_hull = is_it_in_hull(self.points, xyz) 271 in_hull = in_hull.reshape(cube_shape) 272 273 salt_segments = np.zeros(cube_shape, "int16") 274 salt_segments[in_hull == True] = 1.0 275 276 return salt_segments
Insert 3D Salt
Generates a random salt body and insert into input_cube.
generate a random salt body and insert into input_cube salt is modeled using 2 (almost) vertically aligned spheres connected by a convex hull. Points inside the hull are considered salt. center (in x,y) of upper and lower spheres are randomly jittered so salt body is not always predominantly vertical salt body is built using a 'bag of points'. Each point is randomly jittered in x,y,z
Parameters
- top (float): Z-coordinate for top of salt body
- salt_radius (float, optional): Radius of salt body. Defaults to 100.0.
Returns
- None
def
update_depth_maps_with_salt_segments(self) -> None:
278 def update_depth_maps_with_salt_segments(self) -> None: 279 """ 280 Update Depth maps and DepthMapsGaps with salt segments 281 ------------------------------------------------------ 282 283 Updates the depth maps with salt segments 284 285 Parameters 286 ---------- 287 None 288 289 Returns 290 ------- 291 None 292 """ 293 ii, jj = np.meshgrid( 294 range(self.cfg.cube_shape[0]), 295 range(self.cfg.cube_shape[1]), 296 sparse=False, 297 indexing="ij", 298 ) 299 300 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 301 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 302 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 303 304 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 305 306 for ihorizon in range(0, depth_maps.shape[2] - 1): 307 print(" ...inserting salt in horizons...") 308 print( 309 " ...depth_maps min/mean/max = ", 310 depth_maps[:, :, ihorizon].min(), 311 depth_maps[:, :, ihorizon].mean(), 312 depth_maps[:, :, ihorizon].max(), 313 ) 314 depth_map = depth_maps[:, :, ihorizon].copy() 315 depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 316 317 # remove horizon picks inside salt (depth_maps_gaps) 318 depth_maps_gaps_horizon = depth_map_gaps.copy() 319 depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 320 321 _faulted_depth_map = depth_map.copy() 322 faulted_depth_map_indices = _faulted_depth_map.astype("int") 323 faulted_depth_map_indices = np.clip( 324 faulted_depth_map_indices, 325 0, 326 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 327 ) 328 329 _label = salt_segments[ii, jj, faulted_depth_map_indices] 330 depth_maps_gaps_horizon[ 331 depth_maps_gaps_horizon == 0.0 332 ] = np.nan # reset nan's 333 depth_maps_gaps_horizon[_label > 0] = np.nan 334 depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 335 try: 336 number_muted_points = depth_maps_gaps_horizon[_label > 0].shape[0] 337 if number_muted_points > 0: 338 print( 339 " ... horizon {} will have {} points muted inside salt".format( 340 ihorizon, number_muted_points 341 ) 342 ) 343 except: 344 pass 345 346 self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt
Update Depth maps and DepthMapsGaps with salt segments
Updates the depth maps with salt segments
Parameters
- None
Returns
- None
def
update_depth_maps_with_salt_segments_drag(self):
348 def update_depth_maps_with_salt_segments_drag(self): 349 """ 350 Update depth maps with salt segments frag 351 ----------------------------------------- 352 353 Updates depth maps with salt segments and drag. 354 355 Parameters 356 ---------- 357 None 358 359 Returns 360 ------- 361 depth_maps_salt : np.ndarray 362 The depth maps with salt. 363 depth_maps_gaps_salt : np.ndarray 364 The depth maps gaps with salt. 365 """ 366 from scipy import ndimage 367 368 ii, jj = np.meshgrid( 369 range(self.cfg.cube_shape[0]), 370 range(self.cfg.cube_shape[1]), 371 sparse=False, 372 indexing="ij", 373 ) 374 375 depth_maps = self.cfg.h5file.root.ModelData.faulted_depth_maps[:] 376 depth_maps_gaps = self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] 377 salt_segments = self.cfg.h5file.root.ModelData.salt_segments[:] 378 379 if self.cfg.model_qc_volumes: 380 np.save(f"{self.cfg.work_subfolder}/depth_maps_presalt.npy", depth_maps) 381 np.save( 382 f"{self.cfg.work_subfolder}/depth_maps_gaps_presalt.npy", 383 depth_maps_gaps, 384 ) 385 np.save(f"{self.cfg.work_subfolder}/salt_segments.npy", salt_segments) 386 387 depth_maps_salt = np.zeros_like(depth_maps_gaps) 388 depth_maps_gaps_salt = np.zeros_like(depth_maps_gaps) 389 390 relative_salt_depth = 0 391 392 for ihorizon in range(0, depth_maps.shape[2]): 393 print(" ...inserting salt in horizons...") 394 print( 395 " ...depth_maps min/mean/max = ", 396 depth_maps[:, :, ihorizon].min(), 397 depth_maps[:, :, ihorizon].mean(), 398 depth_maps[:, :, ihorizon].max(), 399 ) 400 depth_map = depth_maps[:, :, ihorizon].copy() 401 # depth_map_gaps = depth_maps_gaps[:, :, ihorizon].copy() 402 403 # remove horizon picks inside salt (depth_maps_gaps) 404 # depth_maps_gaps_horizon = depth_map_gaps.copy() 405 # depth_maps_gaps_horizon[np.isnan(depth_map_gaps)] = 0.0 406 407 _faulted_depth_map = depth_map.copy() 408 faulted_depth_map_indices = _faulted_depth_map.astype("int") 409 faulted_depth_map_indices = np.clip( 410 faulted_depth_map_indices, 411 0, 412 self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 413 ) 414 415 _label = salt_segments[ii, jj, faulted_depth_map_indices] 416 if np.max(_label) > 0: 417 relative_salt_depth += 1 418 419 # depth_maps_gaps_horizon[ 420 # depth_maps_gaps_horizon == 0.0 421 # ] = np.nan # reset nan's 422 423 # Apply vertical shift to the horizons where they are inside the salt body, relative to top of salt 424 # depth_maps_gaps_horizon[_label > 0] -= 2 * relative_salt_depth 425 # depth_maps_gaps_horizon = ndimage.gaussian_filter(depth_maps_gaps_horizon, 3) 426 427 # Do the same shift to the non-gapped horizon to build the facies volume 428 depth_map[_label > 0] -= 2 * relative_salt_depth 429 depth_map = ndimage.gaussian_filter(depth_map, 3) 430 431 # depth_maps_gaps_horizon[_label > 0] = np.nan 432 depth_maps_salt[..., ihorizon] = depth_map 433 434 # Re-apply blanking inside the salt to the gapped horizons 435 # _faulted_depth_map = depth_map.copy() 436 # faulted_depth_map_indices = _faulted_depth_map.astype("int") 437 # faulted_depth_map_indices = np.clip( 438 # faulted_depth_map_indices, 439 # 0, 440 # self.cfg.h5file.root.ModelData.faulted_depth.shape[2] - 1, 441 # ) 442 # _label = salt_segments[ii, jj, faulted_depth_map_indices] 443 444 # depth_maps_gaps_horizon[_label > 0] = np.nan 445 # depth_maps_gaps_salt[..., ihorizon] = depth_maps_gaps_horizon 446 447 try: 448 number_muted_points = depth_maps[_label > 0].shape[0] 449 if number_muted_points > 0: 450 print( 451 f" ... horizon {ihorizon} will have {number_muted_points} points muted inside salt" 452 ) 453 except: 454 pass 455 456 if self.cfg.model_qc_volumes: 457 np.save( 458 f"{self.cfg.work_subfolder}/depth_maps_salt_prepushdown.npy", 459 depth_maps_salt, 460 ) 461 462 depth_maps_salt = push_down_remove_negative_thickness(depth_maps_salt) 463 464 for ihorizon in range(0, depth_maps.shape[2] - 1): 465 depth_map_gaps_salt = depth_maps_salt[:, :, ihorizon].copy() 466 faulted_depth_map_indices = depth_map_gaps_salt.astype("int") 467 faulted_depth_map_indices = np.clip( 468 faulted_depth_map_indices, 469 0, 470 salt_segments.shape[2] - 1, 471 ) 472 _label = salt_segments[ii, jj, faulted_depth_map_indices] 473 depth_map_gaps_salt[_label > 0] = np.nan 474 depth_maps_gaps_salt[..., ihorizon] = depth_map_gaps_salt 475 476 depth_maps_gaps_salt[np.isnan(depth_maps_gaps)] = np.nan 477 478 # self.cfg.h5file.root.ModelData.faulted_depth_maps[:] = depth_maps 479 # self.cfg.h5file.root.ModelData.faulted_depth_maps_gaps[:] = depth_maps_gaps_salt 480 481 if self.cfg.model_qc_volumes: 482 np.save(f"{self.cfg.work_subfolder}/depth_maps_salt.npy", depth_maps_salt) 483 np.save( 484 f"{self.cfg.work_subfolder}/depth_maps_gaps_salt.npy", 485 depth_maps_gaps_salt, 486 ) 487 np.save(f"{self.cfg.work_subfolder}/facies_label.npy", _label) 488 489 return depth_maps_salt, depth_maps_gaps_salt
Update depth maps with salt segments frag
Updates depth maps with salt segments and drag.
Parameters
- None
Returns
- depth_maps_salt (np.ndarray): The depth maps with salt.
- depth_maps_gaps_salt (np.ndarray): The depth maps gaps with salt.
def
push_down_remove_negative_thickness(depth_maps: numpy.ndarray) -> numpy.ndarray:
492def push_down_remove_negative_thickness(depth_maps: np.ndarray) -> np.ndarray: 493 """ 494 Remove negative thicknesses 495 --------------------------- 496 497 Removes negative thicknesses. 498 499 Parameters 500 ---------- 501 depth_maps : np.ndarray 502 The depth maps. 503 504 Returns 505 ------- 506 depth_maps : np.ndarray 507 The depth maps with the negative thicknesses removed. 508 """ 509 for i in range(depth_maps.shape[-1] - 1, 1, -1): 510 layer_thickness = depth_maps[..., i] - depth_maps[..., i - 1] 511 if np.min(layer_thickness) < 0: 512 np.clip(layer_thickness, 0, a_max=None, out=layer_thickness) 513 depth_maps[..., i - 1] = depth_maps[..., i] - layer_thickness 514 515 return depth_maps
Remove negative thicknesses
Removes negative thicknesses.
Parameters
- depth_maps (np.ndarray): The depth maps.
Returns
- depth_maps (np.ndarray): The depth maps with the negative thicknesses removed.