GemPy_CrossSection.py

#!/usr/bin/env python
# coding: utf-8

# In[47]:


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


# In[48]:


wells = pd.read_csv(r"C:\Projects\pythonscripts\welldata.csv")  # Change path/file name to match where your CSV actually lives
    #Example Data - ID, well_name, x, y, ground_elevation, formation_name, rock_type, top_depth, bottom_depth, description


# In[49]:


print(wells.head())
print()
print(wells.columns.tolist())
print()
print(f"Number of rows: {len(wells)}")
print(f"Number of wells: {wells['well_name'].nunique()}")
print(f"Formations: {wells['formation_name'].unique()}")


# In[50]:


wells["top_elevation"] = wells["ground_elevation"] - wells["top_depth"]
wells["bottom_elevation"] = wells["ground_elevation"] - wells["bottom_depth"]


# In[53]:


print("Rock types:")
print(wells["rock_type"].unique())
print("\nNumber of intervals by rock type:")
print(wells["rock_type"].value_counts())


# In[54]:


units = wells["rock_type"].unique()

print("Geological units:")
for i, unit in enumerate(units):
    print(i, unit)


# In[55]:


unit_map = {unit: i for i, unit in enumerate(units)}

wells["unit_id"] = wells["rock_type"].map(unit_map)

print(wells[["rock_type", "unit_id"]].drop_duplicates().sort_values("unit_id"))


# In[56]:


# Sort the data by well and depth
wells = wells.sort_values(["well_name", "top_depth"]).reset_index(drop=True)

# Identify where the rock type changes within each well
wells["rock_type_change"] = (
    wells["rock_type"] != wells.groupby("well_name")["rock_type"].shift()
)

# Show only the rows where a new rock type begins
contacts = wells[wells["rock_type_change"]].copy()

print(contacts[
    [
        "well_name",
        "top_depth",
        "bottom_depth",
        "rock_type",
        "top_elevation"
    ]
].to_string(index=False))


# In[58]:


# Keep only rows where the rock type changes
# and exclude the first logged interval in each well

contacts = wells[
    wells["rock_type_change"] &
    (wells.groupby("well_name").cumcount() > 0)
].copy()

# The top of the new interval is the contact elevation
contacts["contact_elevation"] = contacts["top_elevation"]

# Keep only the information we need
contacts = contacts[
    [
        "well_name",
        "x",
        "y",
        "contact_elevation",
        "rock_type"
    ]
].reset_index(drop=True)

print(contacts.to_string(index=False))


# In[59]:


# Start with the complete well data
contact_data = wells.copy()

# Identify the rock type immediately above each interval
contact_data["rock_type_above"] = (
    contact_data.groupby("well_name")["rock_type"].shift()
)

# A contact occurs when the rock type changes
contact_data["is_contact"] = (
    contact_data["rock_type"] != contact_data["rock_type_above"]
)

# Exclude the first interval of each well
contact_data = contact_data[
    contact_data["rock_type_above"].notna() &
    contact_data["is_contact"]
].copy()

# The top of the new interval is the contact elevation
contact_data["contact_elevation"] = contact_data["top_elevation"]

# Keep the useful information
contacts = contact_data[
    [
        "well_name",
        "x",
        "y",
        "contact_elevation",
        "rock_type_above",
        "rock_type"
    ]
].copy()

# Rename the current rock type to make the meaning explicit
contacts = contacts.rename(
    columns={"rock_type": "rock_type_below"}
)

print(contacts.to_string(index=False))


# In[60]:


# Create a clean table of geological contacts
contact_points = contacts[
    [
        "x",
        "y",
        "contact_elevation",
        "rock_type_below"
    ]
].copy()

# Rename columns for clarity
contact_points = contact_points.rename(
    columns={
        "contact_elevation": "z",
        "rock_type_below": "rock_type"
    }
)

print(contact_points.head(20).to_string(index=False))


# In[61]:


print(
    contact_points["rock_type"].value_counts()
)


# In[62]:


import gempy

print(gempy.__version__)


# In[63]:


import gempy as gp
import gempy_viewer as gpv


# In[64]:


surface_points = contact_points.copy()

surface_points = surface_points.rename(
    columns={
        "rock_type": "group",
        "z": "Z"
    }
)

surface_points.head()


# In[70]:


import matplotlib.patches as mpatches

fig, ax = plt.subplots(figsize=(10, 8))

# Colors for each rock type
rock_colors = {
    "Bedrock": "grey",
    "Metamorphic": "darkblue",
    "Sand": "gold",
    "Clay": "brown",
    "Gravel": "orange",
    "Silt": "lightsteelblue"
}

# Give each well a position along the x-axis
well_names = wells["well_name"].unique()
well_positions = {well: i for i, well in enumerate(well_names)}

# Width of each well log
bar_width = 0.25

# Plot each logged interval
for well_name, well_data in wells.groupby("well_name"):

    x_pos = well_positions[well_name]

    for _, row in well_data.iterrows():

        top = row["top_elevation"]
        bottom = row["bottom_elevation"]

        color = rock_colors.get(row["rock_type"], "black")

        rect = plt.Rectangle(
            (x_pos - bar_width / 2, bottom),
            bar_width,
            top - bottom,
            facecolor=color,
            edgecolor="black",
            linewidth=0.5
        )

        ax.add_patch(rect)

# Set X axis
ax.set_xlim(-0.5, len(well_names) - 0.5)
ax.set_xticks(range(len(well_names)))
ax.set_xticklabels(well_names, rotation=45, ha="right")

# IMPORTANT: explicitly set Y limits from the data
y_min = wells["bottom_elevation"].min()
y_max = wells["top_elevation"].max()

padding = (y_max - y_min) * 0.05

ax.set_ylim(
    y_min - padding,
    y_max + padding
)

ax.set_ylabel("Elevation (ft)")
ax.set_xlabel("Well")
ax.set_title("Well Lithology")

# Legend
legend_patches = [
    mpatches.Patch(color=color, label=rock_type)
    for rock_type, color in rock_colors.items()
    if rock_type in wells["rock_type"].unique()
]

ax.legend(
    handles=legend_patches,
    title="Legend",
    title_fontsize=14,
    fontsize=12,
    loc="center left",
    bbox_to_anchor=(1.02, 0.5),
    handlelength=2.0,
    handleheight=1.2,
    handletextpad=1.0,
    labelspacing=1.0,
    borderpad=1.0
)

plt.tight_layout()
plt.show()


# In[ ]: