Skip to content

Create Pseudo Wells ​

This workflow is designed to create pseudo wells within a specified polygon in a Petrel project. The workflow allows users to select a polygon, define the number of pseudo wells, set the TVD (True Vertical Depth) at the first and last points, and choose a folder to place the pseudo wells in. The script generates wells with an adjustable grid approach, ensuring they fit within the polygon and are as evenly spaced as possible. The wells are then created in Petrel with specified wellhead coordinates and surveys.

To learn more about how to create wells, check out this Python Tool Pro tutorial.

Generate Petrel UI ​

python
# Start: PWR Description

from cegalprizm.pycoderunner import WorkflowDescription,DomainObjectsEnum,MeasurementNamesEnum,TemplateNamesEnum

pwr_description = WorkflowDescription(name="Create Pseudo Wells",
                                      category="Wells",
                                      description="Use this workflow to create pseudo wells within an input polygon",
                                      authors="author@company",
                                      version="1.0")

# Use the variable pwr_description to define the UI in the Prizm Workflow Runner and let the Petrel user select the input data.
# This creates a Python dictionary 'parameters' with the GUID and/or values of the user's input data.


pwr_description.add_object_ref_parameter(name='polygon_id',label='Polygon',description='Select a polygon. The polygon must be closed and must contain only 1 polyline,',object_type=DomainObjectsEnum.PolylineSet)
pwr_description.add_integer_parameter(name='pseudo_id',label='Number of Pseudo Wells',description='Define the number of pseudo wells you want to create. ',default_value=200, minimum_value=1, maximum_value=1000)
pwr_description.add_integer_parameter(name='tvd_start_id',label='TVD at first point',description='Define the starting point of the pseudo wells. ',default_value=0, minimum_value=0, maximum_value=10000)
pwr_description.add_integer_parameter(name='tvd_bottom_id',label='TVD at last point',description='Define the last TVD point of the pseudo wells. ',default_value=4000, minimum_value=25, maximum_value=10000)
pwr_description.add_object_ref_parameter(name='folder_id',label='Well Folder',description='Select the well folder you want to place the pseudo wells in.',object_type=DomainObjectsEnum.WellsFolder)


# End: PWR Description

Connect to Petrel and retrieve the user input ​

python
from cegalprizm.pythontool import *
import numpy as np

# Connect to Petrel
petrel=PetrelConnection(allow_experimental=True)
print('PetrelConnection established')

# Retrieve the Polygon selected by the user
polygon = petrel.get_petrelobjects_by_guids([parameters['polygon_id']])[0]

#Check if the polygon is closed and has only one polyline
if len(list(polygon.polylines)) > 1 : 
    raise ValueError(f"{pwr_description.get_label('polygon_id')} : Polygon has more than 1 polyline")

if polygon.is_closed(0) == False:
    raise ValueError(f"{pwr_description.get_label('polygon_id')} : Polygon is not closed.")

print(f"{pwr_description.get_label('polygon_id')} retrieved")

# Assign the user defined number of pseudo wells to a variable
num_wells=parameters['pseudo_id']

# Assign the user defined TVD at first point to a variable
TVD_start=parameters['tvd_start_id']

# Assign the user defined TVD at last point to a variable
TVD_end=parameters['tvd_bottom_id']

# Ensure the start TVD is not larger than the bottom TVD
if TVD_start > TVD_end :
    raise ValueError("Starting depth can't be larger than the TVD at the last point")

# Assign the user selected well folder to a variable
selected_folder=petrel.get_petrelobjects_by_guids([parameters['folder_id']])[0]

#Verify that the well folder has been selected
if selected_folder is None:
    raise ValueError("No well folder has been selected")
print('Well folder retrieved')

Parsing the polygon coordinates ​

The parse_polygon_coordinates function is designed to process polygon data represented as separate lists of x, y, and z coordinates. The function takes these lists (packed in a tuple) as input , unpacks it and returns a list of vertex tuples, which can be used to represent the polygon in subsequent operations, such as plotting or calculating areas.

python
def parse_polygon_coordinates(coord_lists):

    # The input 'coord_lists' is a tuple of three lists: x coordinates, y coordinates, and z coordinates.
    # This line unpacks the tuple into three separate lists and ignores the z coordinates using the underscore '_'.
    x_coords, y_coords, _ = coord_lists

    # This line combines the x and y coordinates into pairs (tuples) representing each vertex of the polygon.
    # 'zip' is used to pair corresponding elements from the x and y lists.
    # The resulting pairs are converted into a list, which represents the vertices of the polygon.
    vertices = list(zip(x_coords, y_coords))
    return (vertices)

polygon_data = polygon.get_positions(0)
polygon_coordinates = parse_polygon_coordinates(polygon_data)
print(f"Polygon coordinates : {polygon_coordinates}")

Approximate the area of the polygon ​

This function is designed to estimate the area of a polygon based on its bounding box. The bounding box is the smallest rectangle that completely encloses the polygon. This method provides an approximation of the area, especially useful when a quick estimation is needed, and the polygon has an irregular shape.

python
def approximate_polygon_area(polygon):
    # Unpacks the polygon's vertices into separate lists of x and y coordinates.
    # The '*' operator is used to unpack the list of tuples into separate tuples for x and y coordinates.
    # 'zip' then pairs each x coordinate with its corresponding y coordinate.
    x_coords, y_coords = zip(*polygon)

    # Calculates the approximate area of the polygon. This is done by multiplying the width and height of the bounding box
    # The width is the difference between the maximum and minimum x coordinates.
    # The height is the difference between the maximum and minimum y coordinates.
    return (max(x_coords) - min(x_coords)) * (max(y_coords) - min(y_coords))

Check if a point is inside the polygon ​

This function uses the ray-casting algorithm to determine if a given point is inside a polygon. The basic idea is to draw a horizontal line from the point to the right and count how many times this line intersects with the polygon's edges. If the number of intersections is odd, the point is inside; if even, it's outside.

Here's what happens in the function:

  • Iterate Through Polygon Edges: The function walks through each edge of the polygon. An edge is defined by two consecutive vertices of the polygon.

  • Horizontal Line Intersection: For each edge, the function checks if a horizontal line from the point would intersect with it. This is done by comparing the y-coordinate of the point with the y-coordinates of the edge's vertices and checking if the x-coordinate of the point is to the left of the edge.

  • Count Intersections: Each time an intersection is found (meaning the point is to the left of the edge and aligned vertically with it), the function toggles the inside variable. This effectively counts how many times the point crosses the boundary of the polygon when moving horizontally to the right.

  • Determine Inside or Outside: At the end of the iteration, if inside is True, the point is inside the polygon; if False, it's outside.

python
def is_point_inside_polygon(x, y, polygon):
    # Get the number of vertices in the polygon.
    n = len(polygon)

    # Initialize a variable to track whether the point is inside the polygon.
    inside = False

    # Get the coordinates of the first vertex of the polygon.
    p1x, p1y = polygon[0]

    # Iterate over the edges of the polygon by walking through each pair of vertices.
    for i in range(n + 1):
        # p2x, p2y are the coordinates of the next vertex in the polygon.
        # The modulo operator (%) is used to loop back to the first vertex at the end.
        p2x, p2y = polygon[i % n]

        # Check if the y-coordinate of the point is between the y-coordinates of the edge.
        if y > min(p1y, p2y):
            # Check if the x-coordinate of the point is to the left of the edge.
            if y <= max(p1y, p2y):
                # Calculate the x-coordinate where the horizontal line through the point intersects the edge.
                # This is only done if the edge is not horizontal.
                if x <= max(p1x, p2x):
                    if p1y != p2y:
                        xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                    # If the point is to the left of this intersection point, it affects the inside/outside status.
                    if p1x == p2x or x <= xints:
                        inside = not inside
                        
        p1x, p1y = p2x, p2y

    return inside

Generate an adjustable grid of points within the polygon ​

This function is designed to place a specified number of wells within a polygon by generating a grid of points and adjusting the grid size to accommodate all the wells.

python
def generate_adjustable_grid(polygon, num_wells):
    
    # Find the minimum and maximum x and y coordinates of the polygon to determine its bounding box.
    minx, miny = min(x for x, y in polygon), min(y for x, y in polygon)
    maxx, maxy = max(x for x, y in polygon), max(y for x, y in polygon)

    # Calculate the approximate area of the polygon using its bounding box.
    area = approximate_polygon_area(polygon)

    # Initial estimate of grid spacing based on the approximate area and the number of wells.
    grid_spacing = np.sqrt(area / num_wells)  

    # Loop to adjust the grid spacing to fit the wells.
    while True:
        # Initialize an empty list to store grid points that fall inside the polygon.
        grid_points = []
        # Create a grid over the polygon's bounding box with the current grid spacing.
        for x in np.arange(minx, maxx, grid_spacing):
            for y in np.arange(miny, maxy, grid_spacing):
                # Check if the current grid point is inside the polygon and if we haven't exceeded the number of wells.
                if is_point_inside_polygon(x, y, polygon) and len(grid_points) < num_wells:
                    # Add the grid point to the list of well locations.
                    grid_points.append((x, y))

        # Check if the number of grid points is sufficient or if the grid spacing is too small.
        if len(grid_points) >= num_wells or grid_spacing < 0.01:  # Minimum spacing threshold to prevent infinite loop
            break
        
        # Reduce the grid spacing to attempt to fit more wells in the next iteration.
        grid_spacing *= 0.95  # Reduce spacing to fit more points

    return grid_points[:num_wells]  # Return only the required number of well locations

# Generating well coordinates on the adjustable grid
adjustable_grid_placed_wells = generate_adjustable_grid(polygon_coordinates, num_wells)


adjustable_grid_placed_wells

Write the pseudo wells back to Petrel ​

python
# Iterate through the list of well coordinates
for wll in range(len(adjustable_grid_placed_wells)):
    # Create new well and place it in the user selected folder
    new_well=petrel.create_well('Pseudowell'+str(wll) ,selected_folder)
    # Set well head coordinates using the generated well coordinates
    new_well.wellhead_coordinates=adjustable_grid_placed_wells[wll]
    # Set the well datum
    new_well.well_datum=('KB',25)
    # Create new deviation survey
    new_surv= new_well.create_well_survey('Pseudo surv','X Y Z survey')
    new_surv.readonly=False
    # Set trajectory of the well. In this case we're creating straight wells (first and last x,y-values are the same). Set the start and end TVD using the user defined values
    new_surv.set(xs=[adjustable_grid_placed_wells[wll][0],adjustable_grid_placed_wells[wll][0]], ys=[adjustable_grid_placed_wells[wll][1],adjustable_grid_placed_wells[wll][1]], zs=[TVD_start,-TVD_end])
    # Set the survey as active
    new_surv.set_survey_as_definitive()
print("Wells successfully created")

create_wells.gif